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abstract 


The  scarcity  and  rising  cost  of  pertroleum  have  motivated 
international  Interest  In  developing  hybrid  automobiles  using  flyuheels 
for  mechanical  energy  storage.  Rlm-type  composite-material  flywheels  are 
prosdslns  designs  for  such  developments.  These  flywheels  significantly 
differ  from  turblne/compressor  systems  in  two  respects.  First,  the  fly¬ 
wheel  rim  attachment  to  Its  hub  Is  very  flexible,  for  both  translation  and 
tUtlng.  Secondly,  these  flexibilities  depend  upon  rotational  speed  through 
centrifugal  stiffening.  In  this  Investigation,  free  whirling,  stability, 
and  forced  wblrllng  are  examined  for  these  flywheel  systems.  The  numeri¬ 
cal  results  presented  here  are  most  directly  applicable  to  Che  Sandla  single- 
rl.  systems  currently  under  development.  However,  the  present  analyses 
can  be  extended  to  other  flywheel  designs  within  the  broad  category  of  the 

rim  type. 

In  the  free-whirling  analysis,  predicted  critical  speeds  are 
encountered  In  the  design  operating  speed  range.  Practical  ways  to  Increase 

such  critical  speeds  are  suggested. 

Effects  of  material  internal  damping  on  the  stability  of  the 

system  are  Incorporated  through  adopting  complex  moduli  in  the  formulation. 

It  is  found  that  the  adverse  effect  of  internal  damping  on  the  onset  ot  in¬ 
stability  can  be  overcome  by  providing  an  adequate  external  damper  up  to 

a  considerably  high  speed. 

Pored  whirling  excited  by  unbalance  and  Initial  tilt  of  rlm 
element  is  studied.  Minimum  external  damping  is  determined  such  that  the 


tnaximum  response  docs  not  exceed  a  permissible  value. 
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1.1.  Historical  Background 

The  scarcity  and  rising  cost  of  petroleum  have  motivated  inter¬ 
national  interest  in  conservation  of  petroleum.  In  view  of  the  large  auto¬ 
motive  use  of  petroleum,  an  important  means  of  accomplishing  this  could  be 
by  introduction  of  hybrid  automobiles  using  flywheels  for  mechanical  energy 
storage.  Most  of  the  flywheels  developed  to  date  have  been  rather  heavy 
and  rigid  ones  constructed  of  materials  that  are  essentially  homogeneous 
and  isotropic,  such  as  steel.  However,  with  the  continuing  development  of 
improved  filamentary  composite  materials,  there  is  considerable  promise  in 
the  use  of  such  advanced  materials  for  flywheel.  The  primary  motivation  for 
this  promise  is  the  high  strength-to-weight  values  for  these  materials, 
since  the  energy  storage  per  unit  weight  can  easily  be  shown  to  be  directly 
proportional  to  this  material  parameter  [l].  Another  advantage  is  the  less 
catastrophic  nature  of  the  failure  of  such  flywheels;  this  can  result  in 
reduction  of  the  weight  of  the  containment  system  which  surrounds  the  fly- 


wheel • 


The  current  national  effort  to  conserve  petroleum  as  an  energy  source 
and  the  current  high  percentage  of  petroleum  being  used  in  automobiles 
motivated  Congress  to  pass  the  Electric  and  Hybrid  Vehicle  Research, ^Develop¬ 
ment  and  Demonstration  Act  of  1976  (P.  L.  94-413).  Subsequently  ERDA* 

*  Energy  Research  and  Development  Administration,  which  was  merged  into  the 
the  Department  of  Energy  on  October  3,  1977. 


established  its  Electric  and  Hybrid  Vehicle  Demonstration  Project  [2].  In 
the  future,  such  a  vehicle  is  a  prime  candidate  for  application  of  composite' 
material  flywheels,  and  the  Department  of  Energy  is  mounting  a  significant 
research  and  development  program  in  this  area,  as  attested  by  many  papers 
presented  at  the  1975  and  1977  Flywheel  Technology  Symposia  [3,  A].  Activ¬ 
ities  of  Sandia  Laboratories  in  this  area  were  reported  as  early  as  1974 
[5],  and  a  major  workshop  on  hybrid  vehicles  was  held  in  1976  [6]. 

Numerous  configurations  of  filamentary-material  flywheels  have  been 
proposed  and  some  have  been  built  and  tested.  Many  of  these  have 
been  assessed  in  a  recent  report  [7]  to  ERDA.  It  suffices  to  mention  here 

the  following  gneral  classes  of  configurations: 

1.  Matrixless  designs  such  as  the  radial  brush  [8]  and  loose- 

fiber-ring-types 

2.  Disk  types,  either  filament-wound  circumferentially  [9] 
or  laminated  [lOJ 

3.  Rim  types,  either  simple  or  multiple,  as  proposed  by  the 
Posts  [ll] 


A.  Shell  types 

The  inherent  flexibility  of  composite-material  flywheels  makes  them 


potentially  vulnerable  to  dynamic  problems  [l2]  which  differ  from  those 
previously  encountered  in  more  rigidly  mounted  metallic  rotors  such  as 
those  in  compressors  and  turbines.  In  this  investigation,  the  dynamic 
phenomenon  of  the  whirling  mode,  recognized  as  one  of  the  most  catastropic 
ones  associated  with  rotating  machinery[ 13-16],  is  carefully  examined. 


To  the  best  of  presently  available  knowledge,  the  only  previous 
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analysis  concerning  the  whirling  dynanics  of  rim-type  flyvheels  existing  in 
the  literature  was  due  to  MciClnnon  [17].  He  considered  a  tvo-mass  system 
having  three  degrees  of  freedom:  hub  tilting,  hub  translation  and  rl.  tilting. 
Ee  did  not  consider  the  radial  flexibility  of  the  connector  with  which  the 
flywheel  rim  is  attached  to  the  hub.  Also,  only  free  whirling  analysis  was 

covered  in  [17]. 

In  Che  present  work,  free  whirling,  forced  whirling,  and  stability 
analyses  are  studied  tor  a  four-mass,  eight  degree-of-treedo.  system  believed 
CO  be  the  most  appropriate  model  to  describe -the  Sandia  flywheel  [l8]. 

Since  the  research  reported  here  was  conducted  in  support  of  an  experi' 
mental  single-thlck-rlng  flywheel  undergoing  research  at  Sandia  Laboratories, 
the  numerical  results  presented  here  are  most  directly  applicable  to  that 
configuration.  However,  with  the  use  of  appropriate  values  of  stiffnesses, 
masses  and  the  dimensions,  many  of  the  analysis  described  would  be  applicable 
to  other  flywheel  designs  within  the  broad  category  of  the  rim  type. 

1.2  nesr-ription  of  the  Sandia  Band-Supported,  Rim^Type  Flywheel 

The  Sandia  Laboratories  experimental  flywheel  prototype  (see  Fig. 
1.1).  for  which  the  analyses  presented  here  are  intended  to  be  applied, 
consists  of  a  thick  rim  of  circumferentially  wound  graphite-epoxy  attached 
to  a  central  aluminum  hub  by  means  of  radially  wound  bands  of  aramid-epoxy. 

The  hub  is  attached  to  the  turbine  through  a  relatively  flexible  steel  fly¬ 
wheel  shaft.  In  turn,  the  turbine  is  supported  by  a  steel  turbine  shaft 


*  It  is  noted  that  some  program  errors  were  fo.und  in  [17] 


TURBINE 


Fig.  1.1.  Schematic  diagram  of  the  Sandia  flywheel  system,  as  installed 
in  the  Sandia-Livermore  spin-test  facility. 
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on  a  set  of  ball  bearings.  An  external  damper  is  located  between  the  lower 
bearing  and  hub  on  the  flywheel  shaft.  The  side  and  plan  views  of  the  fly¬ 
wheel  are  shown  in  Figures  1.2  and  1.3.  It  can  be  seen  that  the  ring  is 
semi-elliptic  in  cross  section  and  the  bands  are  uniformly  rectangular  in 
cross  section* 

The  elastic,  strength,  and  density  properties  of  the  composite 
materials  used  in  the  rim'  and  bands  of  the  Sandia  flywheel  systems  are 
listed  in  Table  1.1. 

In  [18],  Reedy  presented  an  analysis  and  predicted  that  the  design 
energy  storage  goal  of  0.56  kwh  (22.4  watt-hr/lb)  is  achieved  at  a  ro¬ 
tational  speed  of  32,000  rpm  and  that  the  governing  material  strength  is 
reached  at  39,850  rpm. 


1.3. 


In  all  of  the  analyses  presented  here,  the  following  engineering 


assumptions  are  made: 

1.  No  statistical  distributions  of  imperfections  are  considered, 
i.e.  all  eccentricties  are  assumed  to  be  unique  values.  If  sufficient 
statistical  infoinnation  were  to  become  available,  the  approach  recently 
introduced  by  Boyce  and  Kozik  [l9]  could  be  used. 

2.  The  bearings  are  assumed  not  only  to  have  the  same  stiffness 
in  any  diametral  plane,  but  also  to  act  as  simple  supports. 

3.  All  rotating  components  are  assumed  to  be  axisymmetric  in  cross 


section  and  stiffness  properties. 

4.  All  nonlinear  effects  are  neglected  (potential  sources  of  non¬ 
linearity,  include  large  deflections,  bearing  stiffness  nonlinearities. 


bearing  clearances  and  naterial  nonlinearltles. ) 

5.  All  temperature  effects  are  neglected. 

6.  Although  the  flywheel  shaft  is  vertical,  it  is  assumed  that 
the  pendulum  effect  is  small. 

7.  The  shafts  and  bands  ate  relatively  flexible  and  of  relatively 
lev  mass  such  that  they  can  be  modeled  as  discrete,  massless  flexible 
elements. 

8.  The  fli-vheel  ritn,  flj-wheel  hub,  and  turbine  disk  are  assumed  to 
be  rigid  so  that  they  can  be  considered  to  be  discrete  rigid  masses. 

9.  The  svstem  rotates  at  a  constant  angular  velocity. 


It  was  shown  in  [ZO]  that  the  rim  has  relatively  high  frequencies. 


be  treated  as  a  rigid  one. 


Threfore  it  is  justified  that  the  rim  element  can 


CHAPTER  II 


FREE  WHIRLING  OF  THE  UJTDAVPED  SYSTEM 


2.1  The  Gyroscopic  Action  of  a  Rotor 

When  a  shaft  carrying  a  rigidly  attached  rotor,  overhung  with 

respect  to  its  supports,  is  rotated  at  constant  speed,  the  centrifugal 
force  produced  by  an  eccentricity,  regardless  of  how  small  the  eccentric¬ 
ity  is.  causes  the  disk  to  tilt  with  respect  to  the  axis  of  rotation.  This 
tilting  action,  in  turn,  produces  a  gyroscopic  couple  about  the  diametral 
axis  of  the  disk.  The  sense  of  the  gyroscopic  couple  is  such  that  it 
effectively  stiffens  the  shaft  at  the  speeds  associated  with  forward  pre¬ 
cession  and  reduces  the  effective  stiffness  at  the  speeds  associated  with 
retrograde  (or  backward)  precession.  This  gyroscopic  effect  on  rotor 
whirling  was  first  investigated  by  Stodola  [Zl]  in  1918.  see  also  [ZZ]. 

Let  us  consider  a  rotor  which  is  pivotally  mounted  to  a  shaft 
rotating  at  an  angular  velocity  of  n  as  shown  in  Figure  A.l  The  follow¬ 
ing  relationship  can  be  obtained  based  on  a  system  of  fixed  coordinates 
X,  y,  2  (see  Appendix  A  for  detailed  derivation  of  equations  (Z.l)). 


F  =»  -  m  X 

X 


-my 


m  ^x  a  y 


(2.1a) 

(2.1b) 


(2.1c) 


(2. Id) 


lere  (*)  =  a()/3t;  t  is  time;  x  and  y  are  displacements  of  the  rotor  center 

•espectivelyinthe  directions  of  x  and  y;  and  F^  are  the  forces  applied  on 
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-li¬ 


the  shaft  respectively  in  x  and  y  directions  due  to  centrifugal  action 
of  the  rotor;  m  is  the  mass  of  the  rotor;  1^  and  are  the  mas 

moments  of  inertia  of  the  rotor  about  its  respective  axial  and  diametral 
directions;  4  and  4  are  the  angles  of  rotation  of  the  rotor  in  the 

X  j 

respective  xz  and  yz  planes;  and  are  the  moments  applied  on 

the  shaft  respectively  in  the  xz  and  yz  planes.  It  is  noted  that 
small  deflections  and  angles  were  assumed  in  the  derivation  of  equation 
set  (2.1),  i.e.  higher  order  effects  were  neglected. 


If  complex  notation  is  adopted,  the  following  relations  can  be 


obtained: 


where 


F  ■  -  m  r 


F  =  F  +  i  F 

X  y 


M  =  M^  +  iMy 

r  =  X  +  i  y 
^  ^  ^ 


(2.2a) 

(2.2b) 


(2.3a) 

(2.3b) 

(2.3c) 

(2.3d) 


and  i  «  /T  .  It  is  noted  that  the  gyroscopic  effect  was  brought  into 
the  formulation  through  the  second  terms  on  the  right  hand  side  of 
equation  (2.2b). 


2.2  ■  Formulation  of  the  Equations  of  Motion  for  the  Undamped  Flywheel 


System 


In  the  case  of  the  Sandia  thick-rim  flywheel  (see  Figure  1.1), 
the  flywheel  rim  is  connected  to  the  hub  by  means  of  relatively  flexible 
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bands,  in  contrast  to  the  usual  situation  for  overhung  compressor  disks 
which  customarily  are  rigidly  attached  to  the  shaft.  Furthermore,  the 
flywheel  hub  has  an  appreciable  mass  and  has  a  diametral  mass  moment  of 
inertia  which  is  larger  than  its  axial  mass  moment  of  inertia  “ 

0.218).  This  is  the  opposite  of  the  case  of  the  rim  for  which  * 

1.965;  therefore,  the  inertia  couples  for  the  rim  and  hub  will  act  in 
opposite  directions  in  the  range  of  0.218  <  ui/O  <  1.965.  For  ail  of 

the  above  reasons,  it  is  apparent  that  the  rim  and  hub  have  to  be  considered 
as  two  different  masses  rather  than  a  single  disk.  A  few  analyses  of 
multi-mass  disk-shaft,  have  appeared  in  the  literature;  see  for 
instance,  [23],  [24],  [25]  and  [13]  (pp.  225-230).  Unfortunately,  however, 
these  analyses  are  all  applicable  only  to  the  case  when  the  flexible 
members  connect  to  ground.  In  the  case  of  the  Sandia  flywheel  mounted  in 
the  spin-test  facility,  one  flexible  member  (the  bands)  connects  to  an 

otherwise  free  mass  (the  flywheel  rim) . 

In  the  whirling  analysis,  each  rotor  has  two  generalized  displace¬ 
ments  associated  with  translation  (r)  and  rotation  (<))).  In  the  case  of  the 

Sandia  flywheel  as  depicted  schematically  in  Fig.  1.1,  there  are  eight 

.T 

generalized  displacements  {q^}  -  (r^.  r^, 

Here  and  hereafter  the  subscripts  r,  h,  t,  and  d  refer  to  the  rim,  hub 
turbine,  and  external  damper  respectively.  The  compliance  equations  can 

be  written  in  matrix  form  as  follows: 

{q^}  “  ;  i.  j  “  1.  2,  .  .  .  .,  8  (2.4) 

T 

Here  the  generalized  forces  are  (QJ  =  (F^.  F^,  F^.,  F^,  M^) 
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between  rin  and  hub  as  well  as  the  necessary  symmetry  of  the  array  as 
required  by  Maxwell’s  reciprocal  principle,  some  of  the  compliance  terms 
are  repeated  in  the  array,  which  can  then  be  written  as  follows: 


®11  “l2.  ‘  “33  “12  I  “15  “16  '  *17  *18 


“12  °22  I  “l2  V  I  “25  “66  '  “27  “28 

---  -I - ^ - 1--  - 

*13  *12  ,  *33  *12  ‘  *15  *16  ‘  *17  *18 

'  1  I 

«12  “44  i  *12  *44  1  *25  *66  |  *27  *28 

- J - - 1 - - 1 - 

“l5  *25  I  *15  *25  I  *55  *25  '  *57  *58 

I  1  '  . 


*13  *12 


‘16  “66 


a,,  a.,  .  Ogg  0^5  a^g  j  a^g 


H-  -  - - t- 


*17  *27  I  *17  *27  I  *57  *67  1  77  78 

i  ‘  ' 

*18  *28  I  *18  *28  I  *58  *68  I  *78  “88_ 


(2.5) 


with  o 


11  33 


*33 


“44  ^ 


(2.6) 


Here  is  the  band  stiffness  to  resist  relative  in-plane  translation 

between  rim  and  hub,  while  is  the  band  stiffness  to  resist  relative 

out-of-plane  rotation  between  rim  and  hub.  The  values  of  K^^p  and  K^^p 
increase  not  only  with  the  initial  tension  in  the  bands  but  also  with 
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rotational  speed  due  to  the  centrifugal  stiffening  effect.  Detailed 
derivations  of  and  are  presented  in  Appendices  C  and  D  . 

Derivation  of  as  given  in  equation  (2.5)  is  presented  in  Appendix 


Inverting  equation  (2.4).  one  obtains  the  following  relation: 
{Q^}  »  5  i.  j  -  1.  2,  ....  8  (2.7) 


where  [k^]  is  the  stiffness  matrix  and  is  defined  as  follows: 


(2.S) 


Use  of  equation  (2.2)  in  equation  (2.7).  the  following  relation¬ 


ships  are  obtained: 


-  “r  ^ 


-I  4+iJ 

mr  ’^r  mr  r 


a  a 

-  I  .  +  i  j  K 

oh  ^  h  nih  £1 


- 


-I 

mt  ^t  nit  t 


-  “d  ’'d 


md  d  md  d 


(2.9) 


Assuming  that  the  system  whirls  in  a  normal  mode  with  whirling 
frequency  u  ,  one  seeks  the  following  solutions: 


{q^}  -  {q^} 


j  1*1.2.  .  ... .8 


(2.10) 
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where  is  the  amplitude  of  • 

Substituting  equation  (2.10)  into  equation  (2.9),  one  obtains 

the  following  homogeneous  linear  algebraic  equation  set. 


- 

]}  {q^}  -  0 

(2.11) 

where 

8x8 

matrices  with 

zero 

elements 

except  for 

the  following 

elements: 

«11 

e 

“r  ’ 

^22 

“  ^mr  * 

«33 

“  % 

• 

9 

«44 

- 

^mh  ’ 

^55 

-  m^  ; 

«66 

-  ^«t 

5 

“77 

at 

“d  ’ 

^88 

“  ^md  ’ 

(2.12) 

N22 

m 

J  52  ; 

mr 

^44 

-  J  .  n  ; 

oh 

^66 

»  J  ,52 
mt 

» 

^88 

- 

"md" 

Equation  (2.11)  is  referred  to  as  a  generalized  eigenproblem 
[26,  27].  Such  a  problem  can  be  reduced  to  the  standard-type  eigenproblem 

by  using  the  following  definition. 

r-  1  -  ;  (2.13) 

{pj^}  5  uCq^} 

Using  equation  (2.13)  in  equation  (2.11),  one  obtains  the 
following  equation 

[o]  !  [ill  1  I 

- 1 - I  I  (2.14) 

[B2]  I  [Bilj  |{Pi>  1  [^Pi>, 

where  [O]  and  [l]  are  respectively  8  x  8  zero  and  unitary  matrices, 
and  [Bj^],  [B^]  are  defined  as  follows: 

[B,]  =  [M]’"'  m 

[Bg]  H  [M]'"’  [K] 


(2.14) 


(2.15) 
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Equation  (2.14)  is  one  of  the  standard  eigenproblem.  Thus, 
its  frequency  t»  can  be  solved,  at  each  rotational  speed  ^  ,  by  using 
any  existing  eigenproblem  computer  program. 

2.3  Numerical  Results 

The  numerical  values  of  the  various  mass  and  inertia  parameters 
for  two  different  composite-material  flywheel  systems  are  listed  in  Table 
2.1.  Each  of  the  two  systems  have  the  same  rim,  which  is  constructed  of 
hoop-wound  graphite-epoxy  composite  material  and  is  designed  to  achieve 
an  energy-storage  capacity  of  0.56-kwh  at  31,500  rpm. 

System  A  has  flat-band-type  spokes  of  aramid-epoxy  unidirectional 
composite  material  and  is  shown  in  Fig.  1.2.  There  are  six  complete 
bands,  i.e.  the  bands  are  located  30  degrees  apart  around  the  circumference. 
The  bands  are  wound  with  an  initial  tension  of  360  lb  and  pass  through 
and  are  bonded  to  slots  in  the  hub. 

System  B  differs  from  System  A  primarily  in  the  fact  that 
System  B  has  spokes  which  are  wound  flat  on  to  the  rim  but  undergo  a 
90-degree  twist  to  be  wound  on  to  axially  oriented  pins  at  the  hub.  This 
twist  and  a  lower  initial  winding  tension  in  the  bands  (90  lb  for  System 
B)  results  in  higher-in-plane  compliance  of  the  bands.  However,  the 
tilting  compliance  is  smaller* 

Solving  the  eigenvalue  problem  which  is  represented  by 
equation  (2.14)  one  can  obtain  modal  frequencies  as  a  function 
of  the  rotational  speed.  Plots  of  these  relationships  for 
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Systems  A  and  B  are  shown  in  Figures  2.1  and  2-2.  These  figures  are 
plotted  in  log-log  form  in  order  to  get  them  on  the  paper.  All  eight 
modes,  with  forward  and  retrograde  branches  for  each,  are  shown  for  com¬ 
pleteness.  It  should  be  mentioned  that  the  retrograde  branches  correspond 
to  negative  rotating  speeds,  i.e.,  the  direction  of  rotation  of  the 
whirling  is  a  direction  opposite  tothe  direction  of  the  running  speed. 

The  intersections  of  the  u  vs.  2  curves  with  straight  lines 
of  the  form  «  -  n^^  determine  the  so-called  critical  speeds  which 

correspond  to  the  values  of  running  speed  £2  at  which  dynamic  instability 
may  take  place.  The  value  n(-  Wn)  is  called  the  order  of  the  critical 
speed  and  it  is  usually  either  a  positive  or  negative  integer  or  its 
reciprocal.  Lines  corresponding  to  n  *  +  1  and  n  •  +  2  are  shown  in 
Figures  2.1  and  2.2. 

Critical  speeds  of  positive  orders  are  always  associated  with 
forward  precession,  in  which  the  whirling  phenomenon  takes  place  in  the 
same  direction  as  the  shaft  rotation.  In  contrast,  critical  speeds  of 
negative  orders  are  always  associated  with  retrograde  precession,  in 
which  the  whirling  phenomenon  travels  in  a  direction  opposite  to  the 
direction  of  shaft  rotation. 

There  does  not  appear  to  be  any  unanimity  concerning  which 
orders  of  critical  speeds  are  the  most  critical  ones.  For  example. 

Stodola  [21,  22],  Biezeno  and  Grammel  [l3].  and  Hartog  [U]  emphasize  the 
first  order  (n  -  1)  •  Also  ref.  [l3]  claims  that  backward-precession 
critical  speeds  are  less  dangerous  than  forward-precession. ones.  However. 
Yamada  [28]  observed  experimentally  retrograde  as  well  as  forward  critical 
speeds  of  orders  up  to  8.  Recently  Thomson  et  al.  [29]  reported  on 
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2.1  Mode  Identification  for  System  A:  the  number  denotes  the  mode 

number  and  the  suffixes  F  and  R  denote  the  forward  and  retrograde 
precessional  branches,  respectively. 
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Fig.  2.2  Mode  identification  for  System  B: 

the  suffixes  F  and  R  denote  the 
branches,  respectively. 


the  number  denotes  the  mode  number 
;  forward  and  retrograde  precessional 
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experiments  with  single-mass  flywheel  systems.  They  had  difficulty  de¬ 
tecting  the  lower  whirl  phenomena,  although  they  did  detect  one  first- 
order.  lower  mode  critical  speed.  For  the  upper  mode  they  detected  crit¬ 
ical  speeds  of  retrograde  orders  1.  2.  3.  and  5.  on  one  system  and  1  and 
3  on  the  other,  as  well  as  forward  orders  2.  3.  4.  and  5  on  both 

It  is  noted  that  in  Figure  2,1  there  are  shown  twenty-two  first 
and  second-order  critical  speeds  in  the  rotational  speed  range  up  to 
50.000  rpm:  nine  forward  and  thirteen  retrograde.  According  to  the 
reasoning  of  [l3]  and  [14],  the  most  Important  critical  speeds  would  be 
the  ist-order  forward  ones  at  about  700.  7.600.  and  38.100  rpm.  In  con¬ 
trast.  the  experience  of  [17]  and  [29]  suggest  that  higher  mode  forward 
critical  speeds  of  second-order  are  most  dangerous.  It  is  noted  that  the 

amplitude  ratios  and  v^/r^  increase  rapidly  as  the  running  speed 

approaches  the  fourth-mode,  second-order  forward  critical  speed  (approxi¬ 
mately  24,200  rpm).  This  is  consistent  with  comparison  between  the  present 
results  and  the  experimental  results  in  [3l]  as  discussed  in  Chapter  V.  which 
suggests  that  perhaps  the  second  forward  critical  speed  at  the  fourth  mode 

is  the  most  important  one  for  System  A. 

The  curious  triplo-curvaturc  behavior  (two  infUctlou  points)  of 
the  secoud-oohe  forvard  branch  in  the  proxinity  of  the  fourth-node  forward 
branch  (i.  e.  in  the  vicinity  of  «  -  *5.000  cp.,  0  -  23,000  rpn) 

Is  reminiscent  of  the  curve  veering  phenonenon  discussed  by  Leissa  [30] 
and  was  also  found  by  foKInnon  [17],  who  reported  that  an  experimental 

wheel  failed  in  the  vicinity  of  these  conditions. 

The  major  differences  in  System  B  as  compared  to  System  A  are 
ssu^ller  out-of-plane  baud  conpliance  and  larger  hub  mass  and  axial  monent 
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of  inertia  (m^^.  Again  there  are  twenty-tvo  first-  and  second- 

order  critical  speeds  up  to  50,000  rpm:  nine  forvard  and  thirteen  retro¬ 
grade  as  shown  in  Fig.  2.2.  The  first  forward  critical  speeds  are  of  700, 
7,600.  and  38,100  rpm,  while  the  most  important  one  for  the  second-order 
forward  critical  speeds  is  approximately  24,300  rpm.  It  appears  that  the 
calculated  fourth-mode,  second-order  forward  critical  speed  of  24.200  rpm 
was  not  associated  with  the  excessive  amplitudes,  during  spin  tests  (see 
Chapter  V),  in  the  vicinity  of  29,000-30,000  rpm.  The  most  important 
critical  speed  is  less  certain.  However,  there  is  a  distinct  possibility 
that  it  could  be  a  first-order  forward  one  at  the  second  mode. 

One  of  the  objectives  in  designing  the  flywheel  systems  is 
either  to  eliminate  the  critical  speeds,  if  possible,  in  the  operating 
speed  range  or  to  move  them  beyond  the  maximum  operating  speed.  In  the 
Sandia  System  A  and  B  designs,  the  two  most  dangerous  critical  speeds  are 
the  first  order  forward  one  at  the  second  mode,  second 

order  forward  one  at  the  fourth  mode,  (^£^^24*  view  of  the  n  > 

identical  performance  of  Systems  A  and  B  in  clgures  2.1  and  2.2, 
practically  safficiant  to  prasant  numatlcal  results  tor  only  one  systea 
since  these  results  are  applicable  to  another  syste.  oith  an  error  of  less 

than  0*5%. 

Practical  ways  to  increase  are  (1)  to  Increase  the 

area  moment  of  inertia  of  the  flywheel  shaft  ,  (2)  to  decrease  the 

*  Physically,  only  flexural  rotatory  stiffness  is  needed  since 
(fi  )  is  associated  with  the  flywheel-shaft  tilting  mode.  However,  the 
mos^  JLctical  way  to  increase  the  flexural  rotatory  stiffness  is  to  in¬ 
crease  the  shaft  moment  of  inertia,  which  in  turn  increases  the  flexural 
translational  and  coupling  stiffnesses  of  the  shaft. 
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dia»«ral  «..enr  o£  Inertia,  and  (3)  to  increase  the  angle  of 

the  bands  »ith  respect  to  the  plane  of  the  ri.,  d  .  Effects  of  Increasing 
both  1  and  ;  on  ("„>24 

It  IS  slen  that  (1)  nnaffected,  (2)  1»  i-^^bd  by 

approalnately  35  for  each  one-degree  increase  in  ♦  .  and  (3) 

nnlv  1%  when  I  is  doubled*  Effects  of  decreas- 
increased  by  approximately  only  wnen 

I  /o  ^  .  are  shown  in  Fig. 

ing  and  increasing  on  ^‘^cr^24 

2.4.  The  first  Bo  results  for  Fig.  2.3  given  above  remain  applicable. 

Also,  is  increased  by  approaimacely  18%  when  is  reduced  Co 

80!  of  thi  original.  Critical  speed  remains  unchanged  for  the 

above  too  cases  and  for  the  case  even  uhen  I,  and  ;  are  increased 

s_  •  rn  rr  “y  ^  ta  order  to  increase  critical  speed 

similtaneously  as  shown  in  Fig.  2.5.  la  orner  co 

(fl  .  it  is  necessary  to  increase  the  area  foments  of  inertia  of  Che 

£l>.^hrel  shaft  and  turbine  shaft  as  well  as  the  angle  i  simultaneously 
as  shown  in  Fig.  2.6.  When  the  foregoing  three  parameters  are  doubled 
Simultaneously,  critical  speed  in  increased  from  24.200  rpm  to 

be  seen  in  Fig.  2.6. 


31,800  rpa  as  can 


CRITICAL  SPEED  .  RPM  X  10 


-2 


Fig.  2.3  Effects  of  Increasing  flywheel- shaft  area  moment  of  inertia 
and  band  angle  ♦  on  critical  speeds  ^^cr^24* 
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2.4  Effects  of  decreasing  hub  diametral  mass  moment  of  inertia 


I  ,  and  increasing  band  angle  on  critical  speeds 
mh 


and  (a 


CHAPTER  III 


STABILITY  ANALYSIS 


3.1  Internal  Damping  In  the  Shaft  and  Band£ 

As  in  ordinary  (non-rotating)  vibratory  systems,  internal  fric¬ 
tion  (damping)  is  always  present  in' rotating  shafts  if  the  whirling  speed 
is  different  from  the  shaft  rotational  speed.  This  internal  friction  has 
the  effect  of  resisting  shaft  notion  and  causes  dissipation  of  energy. 

It  has  been  well  known  that  internal  damping  has  the  effect  of  reducing  the 
amplitude  of  a  non-rotating  vibrating  system  and  thus  stabilizes  such  a 
system.  In  rotating  machinery,  however,  if  the  rotational  speed  is  above 
the  lowest  first-order  critical  speed,  a  peculiar  phenomenon  sometimes 
can  be  observed,  i.  e.,  the  amplitude  of  whirling  vibration  gradually 
builds  up  and  finally  leads  to  an  unstable  motion.  To  explain  this  kind 
of  phenomenon,  the  imperfection  of  shaft  elastic  properties  must  be  considered. 

There  have  been  two  approaches  used  in  the  literature  to  describe 
internal  danpins:  the  vUcous  type  [32.  33]  and  the  hystereaia  type  [29,  33], 
The  fotmet  is  incorporated  in  a  similar  »ay  as  in  the  simple  sprlng-mass- 
dashpot  system  where  the  friction  force  is  proportional  to  the  velocity  and 
thus  the  energy  dissipation  is  proportional  to  the  frequency.  This  char¬ 
acteristic,  however,  is  not  typical  of  most  shaft  materials  where  the  dis¬ 
sipation  is  of  the  hysteresis  type  and  nearly  Independent  of  the  frequency. 
Qualitative  measurements  of  internal  damping  were  first  carried  out  by  Kim¬ 
ball  [34].  He  loaded  an  end-supported  shaft  and  observed  an  additional  de¬ 
flection.  besides  the  elastic  deflection,  in  a  direction  normal  to  the  loading 
uhen  the  shaft  was  set  in  rotation.  He  also  observed  that  this  additional 
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deflectiotx  is  independent  of  the  rotational  speed  and  is  linearly  de¬ 
pendent  on  the  elastic  deflection,  provided  that  the  elastic  deflection 
is  not  too  large.  In  other  words,  the  ratio  of  this  additional  deflec¬ 
tion  to  the  elastic  deflection  is  equal  to  a  constant  within  the  limita¬ 
tion  of  linear  theory. 

Experiments  with  tension  and  compression  show  that  even  at  very 
low  stress  levels,  materials  are  not  perfectly  elastic.  This 
behavior  is  reflected  in  the  hysteresis  loop  as  shown  in  Figure  3.1.  The 
elastic  behavior  is  represented  by  the  straight  line  AOA,  while  the  width 
of  the  loop  depends  on  the  stress  level  applied  in  the  experiment.  Under 
steady-state  conditions,  the  loop  is  stabilized  as  an  ellipse  with  the 
enclosed  area  equal  to  the  amount  of  energy  per  unit  volume  dissipated 
per  cycle  due  to  hysteresis. 

stress 


strain 


Fig.  3.1  Hysteresis  loop  of  material  under  tension  and  compression. 

consider  .n  ovorhangins  £lywhool-sha£t  system  xhich  is  de£Ueted 
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by  a  force  as  shown  in  Figure  3.2.  The  system  has  a  rotational 

speed  C  ,  and  a  whirling  speed  «.  In  order  to  be  consistent  with  the 

right-hand  coordinate  system  nyt.  !)  and  »  should  be  positive  for  the 
directions  shown  in  Fig.  3.2(b).  For  the  case  !!  >  w  .  any  specific 
fiber  in  the  shaft  rotates  counterclockwise.  When  the  outside  fiber 
moves  inside  through  rotation,  it  is  subjected  to  a  stress  variation 
from  compression  to  tension  and  thus  the  upper  portion  of  the  stress- 
strain  turve  (curve  AU.1)  as  shown  in  Fig.  3.1  should  be  used.  Using 
the  same  argument  one  can  conclude  that  the  lower  portion  (curve  ALA)  of 
Fig.  3.1  should  be  used  to  describe  the  inside  fiber  when  it  moves  out¬ 
side.  Since  Che  upper  portion  has  a  higher  stress  level  compared 
the  lower  portion,  an  additional  force  f;  as  shown  in  Fig.  3.2(b)  is 
generated.  The  ratio  of  F^  to  F^  is  egual  to  a  constant,  t*  . 

under  the  limitation  of  linear  theory.  Since  f;  lags  F^  by  an  angle 

of  -90°.  the  total  force  F  can  be  written  in  complex  form  as  follows; 


F  -  i  F*  “  (1  “  i 

e  e  s  e 


(3.1) 


Soring  that  F^  -  Kr  (K  is  Che  elastic  stiffness  and  is  pro¬ 
portional  to  the  TounI's  modulus  E,) .  one  can  bring  the  internal  damp¬ 
ing  into  the  formulation  by  replacing  E,  used  in  the  previous  chapter 

with  a  complex  modulus 


n  >  u 


(3.2) 


*  The  constant.  is  usually  called  loss  tangent. 
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Flg.  3.2  (a)  Configuration  of  an  overhung  flywheel-shaft  system; 

(b)  Cross-sectional  view  of  shaft. 
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With  a  similar  argument,  one  obtains  the  following  complex 

* 

modulus  for  H  <  w 

E  -  (1  +  i  Y  )E„  ;  ^ 

8  S  S 

The  flywheel  shaft  shown  in  Fig.  3.2  can  also  be  bent  by  applying 
a  moment  on  it.  For  such  a  case,  one  may  use  the  same  philosophy  to 
reach  the  same  conclusions  as  given  in  equations  (3.2)  and  (3.3) .  It  is 
noted  that  the  results  obtained  above  can  be  applied  to  a  rotating  shaft 
under  any  kind  of  supports* 

In  the  Sandia  flywheel  system,  the  band  material  is  Kevlar  49 
aramid/epoxy.  Since  composite  materials  with  polymeric  constituents 
usually  exhibit  more  internal  damping  than  do  the  common  metallic  materials 
one  would  expect  that  the  bands  exhibit  appreciable  internal  damping. 
Although  the  geometric  configuration  for  the  bands  is  different  than  that 
in  the  shaft,  a  similar  conclusion  can  be  obtained.  Thus  the  following 


ec^uatldns  should  be  used  to 

the  analysis. 

bring 

the  internal  damping  in 

the  bands  into 

A 

"b 

- 

(1-1 

> 

>  0) 

(3.4) 

A 

- 

(1  +  1 

y 

<  0) 

(3.5) 

3.2  Formulation  of  the  Equations  of  Motion  for  the  Damped  System 

In  the  Sandia  flywheel  system,  an  external  damper  is  located 
between  the  lower  bearing  and  hub  on  the  flywheel  shaft-  The  external, 
force  acting  on  the  flywheel  shaft  at  the  position  where  the  external 


damper  located  is 
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(3.6) 


Here  C  is  the  viscous  damping  coefficient, 
d 


Using  the  complex  moduli  in  the  shaft  and  band  materials  as 
given  in  the  previous  section  in  conjunction  with  equation  (3.6),  for  the 
external  damper,  one  can  formulate  the  damped  flywheel  system  in  a  similar 
way  to  that  presented  for  the  undamped  system  in  Chapter  II.  Assuming 
that  the  system  whirls  in  a  normal  mode  with  a  solution  form  of  {q^} 

,  one  obtains  the  following  characteristics  equation  similar  to 

equation  (2.11) 


-  [N]a)  -  [K]}{q^}  =  0 


(3.7) 


where  [h]  and  [n]  are  8x8  matrices  with  zero  elements  except  for  the 
elements  defined  in  equation  (2.12)  and  =  i  U  .  [k]  is  the 

complex  version  of  stiffness  matrix  [k]  given  in  equation  (2.8). 

Using  the  same  transformations  as  given  in  equation  (2.13),  one 

obtains  the  following  equation 


[K]  I  i'^1  [N] 


(3.8) 


3.3  Stability  Criterion 


Perhaps  the  most  popular  stability  criterion  used  in  the 
literature  for  linear  analyses  is  the  Routh-Hurwitz  criterion.  However, 
in  view  of  the  large  numbers  of  degrees  of  freedom  considered  here, 
application  of  this  criterion  by  expansion  of  the  determinant  would  be 
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extreaely  laborious.  Here,  the  Routh  concept  is  used  directly  by 
numerically  determing  the  complex  roots  of  equation  (3.8).  Using  an 
existing  routine!  one  can  compute  the  whirling  frequencies  for  a  given 
rotational  speed  .  In  general,  the  whirling  frequencies  are  complex 
numbers,  i.e.  m  +  Iw^  •  Thus,  the  solution  can  be  written  as 

{q  }  =  e^V  .  The  quantity  is  the  damped  frequency; 

while  the  quantity  -  is  the  growth  factor.  It  is  noted  that  if 

_  >  0  ,  the  amplitude  will  grow  with  time  exponentially  and  the  system 

wUl  become  unstable.  On  the  contrary  if  -  Wj  <  0  ,  the  system  is  in  a 
stable  condition.  The  computing  procedure  for  predicting  the  onset  of 
unstable  motion  is  to  start  with  a  given  suitable  rotational  speed  0 
and  compute  all  the  frequencies  then  increase  a  small  amount  of  and 

repeat  the  computation  until  at  least  one  of  the  growth  factors  becomes 

positive* 

3.4  Numerical  Results 

As  in  the  preceding  chapter,  the  numerical  results  presented 
here  are  applicable  to  both  Systems  A  and  B  since  their  differences  are 
less  than  0.5%.  Values  of  Yg  and  Yj,  are  assumed  to  be  0.005  and  0.0172 
respectively.  Effects  of  the  external  viscous  damping  coefficient 
on  the  onset  of  instability  are  shown  in  Fig.  3.3.  If  =  0  ♦  the 

onset  of  instability  occurs  at  the  lowest  first-order  forward  critical 

,  I-,  «f  r  Oess  than  0.25  Ib-sec/in)  is  required  to 

speed.  Only  a  small  amount  of  (.less 

stabilize  uMt.bU  motion  associated  »ltb  the  first  and  third  modes  for 


*  Detailed  programs  are  present  in  Append: 


GROWTH  FACTOR 


•100 


-200 


-30C 


ROTATIONAL  SPEED  RPM 


Fig.  3.3  Effects  of  external  viscous  damping  coefficient  on  the  onset  of 


Instability. 


£2  >  U  .  When  the  value  of  is  greater  than  approximately  0.2  Ib- 

sec/in,  the  onset  of  instability  jumps  up  to  38,100  rpm.  Increasing  the 

value  of  C  has  very  little  effect  on  Increasing  the  onset  of  instability 
d 

beyond  38,100  rpm.  This  implies  that  a  considerable  amount  of  external 
damping  is  required  in  order  to  stabilize  the  system  associated  with  the 
second  forward  mode.  Alternatively,  this  onset  can  be  brought  up  to  a 
fairly  high  rotational  speed  by  (1)  increasing  the  flyhwheel-shaf t  area 
moment  of  inertia  (Fig.  3.4),  (2)  reducing  the  diametral  mass  moment  of 
inertia  of  the  hub  (Fig.  3.5),  (3)  increasing  the  flywheel-shaft  area 
moment  of  inertia  and  band  angle  ^  simultaneously  (Fig.  3.6),  or  (4) 
increasing  the  flywheel-shaft  and  turbine-shaft  area  moments  of  inertia 
as  well  as  the  angle  I  (Fig.  3.7).  However,  the  most  convenient  and 
practical  ways  to  increase  the  speed  at  which  such  unstable  onset  occurs 


in  the  spih^test  set-^up  are 


the  first  three  or  their  combination. 


ROTATIONAL  SPEED  il ,  RPM  x  10 
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Flg.  3.4  Effects  of  increasing  Ig  “<1  ^  stability. 


*lnh''*mh 'basic 


Fig.  3.5  Effects  of  decreasing  and  increasing  * 


Stability* 


Fig. 
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3.6  Effects  of  increasing  ♦  and  I  linmltaneously  on  the 

s 


stability* 
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Fig.  3.7 


,  ^  „  T  -r  I  .  and  1  simultaneously  tin 

Effects  of  increasing  ?  ,  ♦  ^  » 

critical  speeds  ^  ^^cr^24‘ 


CHAPTER  IV 


FORCED  WHIRLING 

4.1  Formulation 

Although  the  rotors  used  for  high-speed  flywheels  can  be  well 
balanced  in  today's  technology,  it  is  not  possible  to  eliminate  the  initial 
unbalance  completely.  This  initial  unbalance  (eccentricity)  causes  resonance 
in  the  neighborhood  of  first-order  critical  speeds.  In  the  absence  of 
damping  and  at  a  fixed  running  speed,  the  resonance  response  will  build  up 
and  finally  cause  the  failure  of  the  system  no  matter  how  small  the  inital 
unbalance  is.  Another  source  of  excitation  for  first-order  critical  whirling 
is  unavoidable  initial  tilt  of  the  flywheel.  For  reasons  stated  in  the 
preceding  chapter,  no  hysteresis  damping  is  produced  at  the  first-order 
forward  critical  speeds.  In  order  to  reduce  the  response  amplitudes,  there 
must  be  some  external  damper  installed  in  the  flywheel  system.  In  the 
Sandia  flywheel  system,  there  is  an  external  damper  located  on  the  flywheel 
shaft  as  shown  in  Fig.  1.1.  This  external  damper  not  only  effectively 
reduces  the  response  amplitudes,  but  also  is  effective  in  overcoming  adverse 
internal  damping  to  stabilize  the  system  when  the  operating  speed  is  beyond 
the  lowest  first-order  critical  speed.  In  the  present  chapter,  the  frequency 
response  of  the  system  will  be  studied.  Minimum  damping  will  be  determined 

for  a  given  maximuni  response  allowance* 

If  a  rotor  is  unbalanced  with  an  eccentricity  of  e,  the  trans¬ 
lational  displacement  of  the  rotor  gravity  center  r^,  can  be  expressed 

as 
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where  m  is  the  rotor  mass. 


Similarly,  the  tilting  angle  of  a  rotor  with  an  initial  tilt 


can  be  written  as 


•^C  *  •f'  +  ® 


(4.2) 


It  is  noted  that  r  and  ^  in  the  preceding  two  equations  are 
respective  shaft  (or  bands,  if  appropriate)  translation  and  slope. 

A  derivation  similar  to  that  for  equations  (2.2a)  and  (2.2b) 
gives  the  following  two  equations 

,  “  (4.3) 

F 

«.  -.i  -i.-tTnA  (4.4) 


M  '  Vc  ^  ^  •'m''  ' 

where  F  and  M  are  respective  force  and  moment  applied  on  the  shaft. 

Inserting  equations  (4.1)  and  (4.2)  into  equations  (4.3)  and 
(4.4),  respectively,  one  obtains 

.  _  o2>C  (4.5) 


„  .  .1^ ; .  i  n ; .  (i„  -  V  o' . 

The  shaft  at  the  position  where  external  damper  is  located  receives 


the  following  force  : 


“  “d  ’^d  ■  ^d  ’^d 


The  compliance  equations  can  be  written  in  matrix  form  as 


follows ; 


{qf}  “  i.  j  “  1,  2,  .  .  .,  7 

{<l^}  -  *^r*  ^b*  ^d*  ^d^ 

(Qj)  5  {F^.  F^.  F^.  F^-  »d'’ 


(4.8) 


where 
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ani!  are  given  in  Appendix  E 


Due  to  the  geometric  configuration  of  the  Sandia  Flywheel  System 
it  is  expected  that  unbalance  and  initial  tilt  mostly  occur  in  the  rim. 

In  view  of  this,  frequency  response  excited  by  such  kinds  of  geometric 
imperfections  are  studied.  However,  such  studies  could  be  extended  to 
incorporate  similar  excitation  received  from  other  components,  if  necessary 


or  Important. 


Inserting  equations  (4.5)  and  (4.6)  into  equation  (4.8).  and 


using  the  following  solution  form 


{q.}  -  {q^}  e 


one  obtains 


{[Kij]  +  {qi>  =  tfjJ 


(4.9) 


(4.10) 


where  [K^]  S  and  have  zero  elements  except 

the  following 

5  ^22  ~  ^®r  ~  ^mr  ’  *^33  °h 


M44  ”  ^mh  -  -^mh  5  M. 


m  •  M  =*  I  “  d  .  (4.11) 

“c  »  “66  mt  mt  ' 


^7  “  “d  ’  ^88  ’  ’  ^77  *  ^^d» 

{fj}  is  the  excitation  vector  and  is  given  below 
{f^}  -  (e^m^  .  0*  0* 

Equation  (4.10)  is  a  set  of  nonhomogeneous  linear  algebraic 
equations  in  terms  of  the  responses. 
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4.2  Numerical  Results 

Solving  e<,uation(4.10)  on  an  IBM  370,  Model  158J  digital  co.- 
puter,  one  obtains  the  responses  as  a  function  of  the  rotational  speeds. 

The  responses  induced  by  rim  eccentricity  are  presented  in  rigures  4.1  and 
4.2  for  three  different  viscous  damping  coefficients.  For  ■  0  (in 
units  of  Ib-sec/in.  hereafter),  the  amplitudes  become  unbounded-  It  is  inter¬ 
esting  to  see  that  excellent  agreements,  as  far  as  the  first-order  critical 
speeds  concerned,  are  obtained  between  Figures  4.1,  4.2  and  Figure  2.1  of 
Chapter  II.  It  can  be  seen  that  larger  responses  are  obtained  in  the 
lowest  critical  speed  compared  to  that  in  the  two  otherwise  higher  critical 
speeds.  However,  the  former  is  less  important  since  the  design  operating 
speed  is  in  the  range  of  8,000  rpm  and  up.  The  responses  induced  by  rim 
initial  tilt  are  presented  in  Figures  4.3  and  4.4.  Again  good  agreements 
are  obtained  between  Figures  4.3,  4.4  and  Figure  2.1.  In  contrast  to  the 
response  curves  presented  in  Figures  4.1  and  4.2,  the  larger  responses^ 

occur  in  the  third  critical  speed. 

I»  the  Sandle  deslgu,  there  is  .  "stop"  installed  Im  the  system. 

uhlch  limits  the  hub  translation  to  a  maximum  of  0.0283  Inches,  to  order 
to  prevent  the  hub  from  hitting  the  "stop",  an  adequate  amount  of  external 
damping  should  be  present.  According  to  the  Information  provided  by 
c-a..  laboratories,  the  maximum  rim  eccentricity  and  initial  tilt  pertlent 
to  their  designs  are  respectively  0.028  Inches  and  0.15  degrees,  to  other 
vords,  the  maximum  allouable  ratios  of  the  hub  translation  to  rim  eccentricity 
and  initial  tUt  are  respectively  1  and  12  In/rad.  A  plot  of  the  ratio 
of  the  hub  translational  amplitude  to  rim  eccentricity  versus  rotational 

speed  is  shorn,  in  Figure  4.3.  It  can  be  seen  that  the  minimum' amount  of 
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Fig.  4.2  Ratios  of  tilt  amplitude  to  rim  eccentricity  for  »  0,  3 


and  5  Ib-sec/in. 
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Fig.  4.4  Ratios  of  tilt  amplltitudes  to  rim  initial  tilt  for  -  0, 
3,  and  5  Ib-sec/in. 
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Fig.  4.5  Ratio  of  hub  ttanslatioaal  amplitude  to  rim 


eccentricity. 
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viscous  damping  coefficient  required  is  0.25  Ib-sec/in  (as  discussed  in 
the  preceding  chapter,  this  value  is  sufficient  to  maintain  a  stable 
motion  up  to  38,100  rpm) .  As  a  check,  a  plot  of  the  ratio  of  the  hub 
translational  amplitude  to  rim  initial  tilt  is  shown  in  Figure  4.6.  It 
can  be  seen  that  the  peak  responses  for  “  0.25  lb-sec/ in  are  well 
below  the  maximum  allowance  (12  in/ rad) .  The  hub  translational  responses 
induced  by  the  rim  eccentricity  and  initial  tilt  versus  rotational  speed 
are  given  in  Figure  4.7.  Again,  a  value  of  0.25  Ib-sec/in  for  is 

sufficient  to  prevent  the  hub  from  hitting  the  stop  . 

A  conservative  estimation  of  viscous  damping  coefficient  is 

■presented  in  Appendix  F.  Considering  the  worst  case  as  shown  in  Figure  4.7. 

at  n  =  38,100  rpm,  the  corresponding  damper  lateral  translation  is  0.051  in. 

The  damper  clearance  pertinent  to  the  Sandia  spin-test  facility  is  0.040  in. 

This  gives  a  value  of  more  than  unity  for  the  eccentricity  ratio;  which  is 

not  physically  possible,  since  the  eccentricity  ratio  can  not  exceed  a 

value  of  one.  However,  according  to  the  calculation  conducted  in  Appendix  F, 

the  viscous  damping  coefficient  should  exceed  0.25  Ib-sec/in  when  e  »  0.65 

or  a  damper  displacement  of  about  0.026  in.  For  an  eccentricity  e  -  0.9, 

Appendix  F  lists  =  1.3  Ib-sec/in.  Since  this  value  is  considerably 

greater  than  the  0.25  Ib-sec/in  value  on  which  Fig.  4.7  is  based,  it  is  clear 

that  the  0.051-in.  damper  displacement  mentioned  above  would  be  considerably 

reduced  for  C  =1.3  Ib-sec/in.  and  there  would  be  no  danger  of  the 
d 

damper  hitting  its  housing. 
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tilt. 


CHAPTER  V 


COMPARISON  WITH  EXPERIMENTAL  RESULTS 

la  the  epln  tests  of  the  tvo  Sandle  tly»heel  systems  (denoted 
here  as  Systems  A  and  B)  tested  by  Sandla  personnel  at  the  SandlA  Uver- 
note  Lahoratory  [3l].  n  and  y  cootdlnates  of  the  huh  position  la  the 
horizontal  pUne  uere  picked  up  by  prorimity  gages  and  displayed  on  an 


oscilloscope. 

The  System  A  flywheel  displayed  considerable  vlbratioual  amplitude 
in  the  low-speed  range  below  about  2,000  rpn.  but  It  was  possible  to 
accelerate  through  this  speed  range.  (The  excessive  amplitude  may  have 
been  indicative  of  the  predicted  critical  speeds  in  the  vicinity  of  700 
and  2,000  rpm.)  Starting  at  about  14,000  rpm,  there  was  a  gradually 
increasing  buildup  in  amplitude  which  accelerated  rapidly  starting  at 
18,000  rpm.  Expeclally  sudden  increases  in  amplitude  were  noted  at  20,200 
rpm  and  20,900  rpm.  These  sudden  jumps  while  the  rotational  speed  was 
gradually  Increased  appeared  to  be  indicative  of  the  nonlinear  jump  phe¬ 
nomenon  associated  with  a  "softening"  testoring  force  [35].  It  is  also  pos- 
slble  that  these  sud(ien  increases  were  due  to  any  one  of  these  causes.  (1) 
sudden  failure  of  a  band,  (2)  slippage  of  the  bands  in  the  vicinity  of  the  hub, 
or  (3)  rapid  opening  of  a  delamination  in  the  rim.  The  flywheel  shaft  broke 
due  to  excessive  amplitude  when  a  speed  of  approximately  22,100  rpm  was  reached 


By  ultrasonic  and  radiographic  moans  it  had  boon  found  prior  to 
the  spin  tost  that  this  particular  flywheol  rim  appoarod  to  have  some 
locallzod  areas  of  fiber  buckling.  An  earlier  tost  of  another  flywheel  of 
the  System  A  design  had  achieved  17,900  rpm  before  the  lead  balance  weights 
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were  thrown  off  and  the  test  stopped.  Surface  flaws,  but  no  significant 
internal  flaws,  were  detected  in  it  prior  to  test,  and  post-failure  exam¬ 
ination  showed  that  it  had  separated  in  the  vicinity  of  the  surface 
irregularity. 

In  summary,  it  appears  that  the  calculated  fourth-mode,  second- 
order  forward  critical  speed  of  24,200  rpm  for  System  A  is  in  good  agree¬ 
ment  with  the  observed  failure  due  to  excessive  amplitude  at  22,100  rpm. 
This  agreement  was  surprising,  since  all  of  the  stiffness  values  used 
were  calculated  ones,  some  of  which  were  based  on  boundary  condition  assum- 
tions  which  may  not  have  been  sufficiently  realistic. 

The  System  B  flywheel  exhibited  a  generally  smoother  ride  (less 
vibrational  amplitude)  than  did  System  A.  A  possible  explanation  for  this 
could  be  the  inherently  more  s>-mmetric  stiffness  distribution  of  the  band 
geometry  in  System  B.  Pronounced  sudden  decreases  in  amplitude.  Indicative 
of  nonlinear  jump  phenomena  associated  with  a  "hardening"  restoring  force 
[35]  were  observed  at  approximately  22,000  rpm  and  24.800  rpm.  There  was 
a  gradual  increase  in  amplitude  starting  at  about  27.500  rpm,  with  sudden 
increases  in  amplitude  at  approximately  29,000  rpm  and  30,000  rpm.  There 
was  a  loud  report  associated  with  the  first  of  these,  and  final  failure  was 
at  30,100  rpm.  Post-failure  examination  of  this  flywheel  indicated  that 
the  rim  still  retained  its  structural  integrity,  but  many  of  the  bands  had 
failed  near  their  attachement  ot  the  rim.  It  is  not  known  which  of  the 
following  phenomena  directly  caused  failure  of  these  bands:  (1)  snapping 
of  some  of  the  bands  due  to  excessive  steady  centrifugal  and  aerodynamic 
loads,  (2)  excessive  dynamic  mechanical  loads  resulting  from  the  system 
dynamics  causing  the  bands  to  snap,  or  (3)  excessive  abrasion  due  to 
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rubbing  after  the  flywheel  dropped  into  the  spin  pit  following  failure 
of  the  reduced-section  breakaway  shaft.  An  earlier  test  of  the  System  B 
design  had  displayed  dynamic  instability  at  29,000  rpm  and  the  test  was 
discontinued  at  that  point. 

In  summary,  it  appears  that  the  calculated  fourth-mode,  second- 
order  forward  critical  speed  of  24,250  rpm  for  System  B  was  not  associated 
with  the  excessive  amplitudes,  during  spin  tests,  in  the  vicinity  of  29,000- 
30,000  rpm.  However,  there  is  a  calculated  first-order  forward  critical 
speed  for  the  second  mode  (associated  with  flywheel  shaft  rotary  flexibility) 
at  approximately  38,000  rpm.  It  may  be  tentatively  conjectured  that  the 
System  A  flywheel  failed  at  a  second-order  critical  speed  (22,000  rpm)  due 
to  stiffness  asymmetry  of  the  bands,  while  System  B,  could  have  gotten 
beyond  this  speed  only  to  encounter  excessive  amplitude  at  the  next  higher 
first-order  critical  speed  (30,000  r?m) . 


CHAPTER  VI 

CONCLUDING  REMARKS  AND  SUGGESTIONS  FOR  FURTHER  RESEARCH 

In  this  investigation,  free  whirling,  forced  whirling,  and 
stability  of  two  rim-type  composite-material  flywheel  systems  currently 
under  development  and  spin  test  at  Sandia  Laboratories  were  studxed* 

Critical  speeds  were  encountered  in  the  design  operating  speed 
range  (8,000  rpm  to  31,500  rpm) .  For  System  A,  perhaps  the  most  important 
critical  speed  is  the  second-order  forward  one  at  the  fourth  mode  (approx¬ 
imately  24,000  rpm).  Practical  way  to  increase  such  a  critical  speed 
beyond  the  desired  maximum  operating  speed  is  to  increase  the  band  angle 
and  the  flywheel-shaft  and  turbine-shaft  area  moments- of  inertia  (or  to 

decrease  the  lengths  of  the  these  shafts  since  the  purpose  is  to  increase 

/ 

the  shaft  rotatory  stiffness)  simultaneously.  For  System  B,  the  most 
important  critical  speed  is  less  certain.  However,  there  is  a  distinct 
possibility  that  it  could  be  a  first-order  forward  one  at  the  second  mode 
(approximately  38,100  rpm).  This  critical  speed  can  be  moved  up  fay  either 
reducing  the  hub  diametral  mass  moment  of  inertia  or  increasing  the  fly- 
wheel-shaft  area  moment  of  inertia  (or  reducing  the  length  of  this  shaft). 

The  adverse  effect  of  material  internal  damping  on  the  system 
stability  can  be  overcome  by  providing  an  adequate  external  damper.  Only 
a  small  value  of  viscous  damping  coefficient  is  required  to  stabilize  the 
system  up  to  the  first  critical  speed  at  the  second  mode.  However,  a 
considerable  amount  of  external  damping  is  required  to  stabilize  the  sys¬ 
tem  above  this  critical  speed.  Alternatively,  this  onset  of  instability 
can  be  moved  up  by  either  reducing  the  hub  diametral  mass  moment  of 
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Inertia  or  Increasing  the  area  moment  of  inertia  (or  reducing  Che  length) 
of  the  flywheel  shaft. 

Jq  order  to  prevent  the  hub  translational  motion  from  hitting 
the  "stop"  during  forced  whirling  excited  by  unbalance  and  initial  tilt, 
a  minimum  viscous  damping  coefficient  of  0.25  lb~sec/in  is  required  for  the 
Sandia  flywheel  designs. 

Several  new  questions  have  arisen  from  this  research  and  it  is 
recommended  that  Che  following  additional  research  be  conducted: 

1.  In  view  of  the  "nonlinear"  jump  phenomenon  observed  during 
spin  test,  nonlinearity  arising  from  the  band  action,  believed  to  be  the 
most  significant  one  in  the  systems  investigated,  should  be  examined. 

2.  There  exists  a  possiblity  for  the  band  on  the  compressive 
side  of  the  bend  to  buckle  in  the  double-band  action  associated  with  rim 

tilt  action  relative  to  the  hub. 

3.  The  stacking  sequence  of  bands  which  causes  unsymmetric 
stiffness  distribution  circumferentially  may  be  significant  to  the  stability 
of  the  system.  An  investigation  on  this  unsymmetric  band  stiffness  dis¬ 
tribution  may  provide  some  Information  pertinent  to  the  excitation  for  the 
second-order  critical  speeds  in  System  A. 
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APPENDIX  A 


DERIVATION  OF  ROTOR  BASIC  EQUATIONS  INCLUDING  GYROSCOPIC  EFFECT 


Consider  an  axisyinmetrlc  rotor  as  shown  in  Fig.  Al.  The  rotor 
rotates  at  an  rotational  speed  fJ  with  respect  to  axis  which  passes 

through  the  geometric  center  C  of  the  rotor.  Due  to  the  symmetric  nature 
of  the  rotor,  coordinate  system  x^^  y^^  coincides  with  one  of  the  prin- 
cipal  axes  of  the  rotor.  Thus,  the  angular  momentum  H  of  the  rotor 
based  on  system  Xj^  y^^  is 

H  *10  i  H  *1^  ;  (A.~l) 

n  i  _  f  m  y  2  ™ 


J  G  (A-1) 
m 


Here  I  >  mass  moments  of  inertia  of  the  rotor  respectively 

m  m 

about  its  dlameteral  and  axial  directions,  and  are  rotations  of 


the  rotor  respectively  in  the  xz  and  yz  planes,  and 

are  components  of  H  in  the  *1^1  )/3t 

The  angular  momentum  of  the  rotor  based  on  system  x'y'z*  can 
be  obtained  as  follows  by  resolving  the  components  in  equations  (A-1)  and 
using  the  small-angle  approximation. 


(A-2) 
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To  the  first-order  approximation  (small  <t)^  and  the  Euler  s 

equations  of  motion  “  H^y^.  are  applicable  and  hence  the  following 

inertia  couples  are  obtained 
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Equations  (A-3)  represent  the  inertia  couples  applied  on  the 
rotor.  The  inertia  couples  applied  on  the  shaft  by  the  rotor  are  the 
reactions  received  as  of  equations  CA-3).  Thus,  the  nonzero  couples  applied 


on  the  shaft  are 


\  *x  -  ^ 


My  “  "  \  +  V  ^ 


(A-4) 


The  first  terms  on  the  right-hand  sides  of  equations (A-4)  are 
associated  with  rotatory  inertia,  while  the  second  terms  are  associated 
with  gyroscopic  effects. 


appendix  b 


INTERACTION  BETWEEN  THE  RIM  AND  BANDS 

In  this  analysis,  only  the  additional  forces  due  to  centrifugal 
action  are  determined.  They  add  to  the  initial  forces  due  to  prestress 
during  winding. 

The  free  centrifugal  expansion  of  a  single  band  is  calcualted  as 
follows  (see  Fig.  B1  for  detailed  geometry).  First,  the  centrifugal  force 
developed  at  any  arbitrary  radius  r^^  is  determined  by 

p  .  Ajj  0^  r  dr/cos  i  -  4  P^, 


where  R^  is  the  rim  inside  radius;  4  is  defined  in  Fig.  Bl. 
Thus,  the  axial  strain  of  the  band  at  rj^  is 

e^(ri)  *  Pjj^  cos(p/A^E^  =  (Pj^  C  /2  E^^)  (R^  - 


(B-2) 


The  axial  displacement  at  the  outer  edge  of  the  band  (i.e.  at 
its  point  of  attachment  to  the  rim)  is 

.  .  e.(rp<lrj  -  V® 

°°  Jr, 

“  (B-3) 


where  R^  is  the  hub  outside  diameter. 

For  the  rim,  the  hoop  tension  due  to  free  centrifugal  action  is 


p  A  R^ 


(B-4) 


;us,  the  radial  displacement  of  the  rim  at  its  centroidal  radius  R  is 
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-66* 


2  3 

U  -  p  a  R  /E 
rc  r 


(B-5) 


Comparison ‘of  the  values  of  cos  ^  (this  is  the  band  radial 

displacement  at  the  outer  edge)  and  using  data  pertinent  to  the 

present  investigation,  shows  that  the  latter  is  greater  than  the  former. 
Thus,  the  bands  are  loaded  in  tension.  To  determine  the  interaction 
between  the  two,  we  consider  a  typical  "repeating  section"  as  shown  in 


Fig.  B2, 


A  summation  of  radial  forces  yields  the  following  expression 


(Rfj^  -  N“)  sin  6^  »  Q  cos  0^ 


(B-6) 


where  =  circumferential  force,  Q  =  shear  force,  and  f^^  is  the 
centrifugal  body  force  per  unit  length  given  by 


(B-7) 


The  elastic  boundary  conditions  at  the  edges  of  the  segment 


can  be  written  as 


N®  “  K^u, 


^  ^  2  hvK  ~  V 


(B-8) 


where  u  is  the  radial  displacement  of  the  edge  at  8  =  0^  and  the 

o 


stiffnesses  are  given  by 


(B-9) 


where  is  the  unsupported  spaa  length  of  the  band. 

Combining  equations  (3--3)*  (B-6)  through  (B-9),  one  obtains 


the  following  expression 

2  pAR^  tan  ^  ^b  ^  cos^  5  [(R^/ 3 -( R^/2)+(R^/6)]i^2L.^ 

u  a  - - - — - - -  ■  ‘ 

o  2  “ 

(AE/R)  tan  d  -5“ 


(B-10) 
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The  additional  radial  force  in  a  band  due  to  the  rim-band 


interaction  is 


■  T  -  “bo  ““ 


(B-11) 


The  additional  tension  in  a  band  due  to  the  rim-band  interaction 


P  *  «  2  Q  cos  if  =  *S)r^%  ”  “bo  <')cos  f 


(B-12) 


The  total  tension  in  a  band  is  the  sum  of  cos  f  and 

p  '  .  The  former  varies  from  zero  at  the  outer  edge  to  a  maximum  value 
be  . 

at  the  inner  edge.  The  integrated  average  of  P  cos  6  is 


(B-13) 


Performing  the  integration,  one  can  rewrite  equation  (B-13)  as 


S)^  (2  rJ  -  3  eJ  -  V 


(B-14) 


Thus,  the  total  centrifugal  tension  in  a  band  is 


^  +  -be’ 


(B-15) 


APPENDIX  C 


DETERMINATION  OF  IN-PLANE  BAND  STIFFNESS 


A  small  in-plane  displacement  (see  Fig.  Cl)  of  the  rim 
relative  to  the  hub  causes  each  band  set  (that  is,  the  upper  and  lower 
bands  which  are  shown  superimposed  in  the  figure)  not  only  to  stretch 
along  its  axial  direction  but  also  to  bend  in  the  plane  containing  the 
displacement  r^  .  The  former  produces  the  in-plane  stretching  rigidity; 
while  the  latter  produces  the  in-plane  flexural  rigidity. 


The  amount  stretched  in  the  band  set  is 

u.  »  r  cos  (ij).  +9)  cos  0 
j  c  J 


(C-1) 


where  0  is  the  angle  between  the  displacement  direction  and  the  nearest 
band  having  additional  tension  and  is  the  angle  (measured  in  the 

same  direction)  between  that  nearest  band  and  the  j  band. 

Thus  the  axial  force  in  a  band  set  is 


P  =  2(P,  +  P  +  u.) 

Aj  '  i  c  j 


A^  E 

p  =  ZfP.  +  P  + —  r  cos(ii).  +  9)  cos  <()] 
Aj  •-  i  c  L,  c  J 


(C-2) 


where  both  P.  and  P  are  assumed  to  have  the  same  value  in  all 


directions. 


The  force  component  in  the  direction  of  the  displacement  r^  is 

P  =  p  cos  Oil.  +  6)  cos  (}i 
Hj  Aj  J 

-  2  [P  cos  (it;  +  9)  +  P  cos(t(;.  +  9)  +  7—  ^  cos  (ij;  +  9)  cos  ^.Jcos 

*-i  J  cj  S  ^ 

(C-3) 
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CENTER  OF 
ROTATION 
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The  resultant  restoring  force  is 


— r  cos^  0  E  cos  (yi  9) 


p  «  r  P  ■  ■  2  — —  r  cos  (ji  i 

j-1 


A.  E  2  ” 

p  =12  r  cos  ^ 

®  L 


(C-4) 


Thus,  the  Integrated  stiffness  due  to  the  variation  of  axial 


tension  in  the  bands  is 


12  tp.  c=s^ ; 
h) 


fC-5) 


The  in-plane  flexural  rigidity  in  the  bands  is  the  ratio 

of  lateral  force,  applied  on  one  end  of  the  bands,  to  the  side-sway  de¬ 
flection.  Using  classical  beam-tension  bar  theory,  one  obtains 

-1 

=  (kP  sinh  kL^)C2(l  -  cosh  kL^>  +  kL^  sinh  kl^]  (C-6) 

where  P  is  the  total  tension  in  each  band;  k  s  (P^Vl^  *  ^1 

the  area  moment  of  inertia  of  band  cross-sectional  area  about  its  cen- 

troidal  axis  parallel  to  the  flywheel  shaft. 

A  set  of  bands  parallel  to  the  direction  of  load  application 
can  not  carry  any  flexural  load  and  thus  does  not  contribute  to  the 
integrated  flexural  stiffness.  However,  a  set  of  bands  at  90“  to  the 
direction  of  loading  contributes  its  in-plane  flexural  rigidity.  For 
any  given  loading  direction,  the  integrated  in-plane  flexural  rigidity 

is  given  by 
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tegrated  ”  ^Ccos^e  +  cos^  (30"  -  0)  +  cos^(30"  +  0) 

+  cos^(60"  -  0)  +  cos^(60"  +  0)  +  cos^  (90"  -  0)] 

«  T>  V  (C-7) 


The  total  Integrated  in-plane  stiffness  is 

hip  ■ 


(C-8) 


In  the  Sandla  System  A  flywheel,  the  average  is  approximately 
equal  to  2.66".  Using  equation  (C-8),  one  obtains  only  0.2S  less  stiffness 
by  incorporating  the  angle  $  as  compared  to  the  case  when  it  is  not 
incorporated. 

For  System  B,  the  bands  are  wound  flat  on  to  the  rim  but  undergo 
a  90-degree  twist  to  be  wound  on  to  axially  oriented  pins  at  the  hub. 

Bands  in  such  systems  are  shewn  in  Fig.  C2(a).  Application  of  equation 
(C-8)  is  not  valid  unless  a  modification  is  carried  out.  Physically  the 
first  term  on  RES  of  equation  (C-8) ,  which  is  associated  with  stretching 
in  axial  direction,  is  unaffected.  However,  the  second  term  should  be 
decreased  due  to  reduction  of  moment  of  inertia  along  x  . 

If  the  force  applied  at  the  right-hand  side  is  Q  as  shown  in 
Fig.  C2(a),  the  moment  distribution  is 

M(x)  -  -  Q(Ljj  -  x)  (C-9) 


Using  the  same  philosophy  as  presented  in  [36],  the  trans¬ 
lational  deflection  at  the  point  where  Q  is  applied  can  be  expressed  as 
follows  by  noglficting  lbs  she^r  floxibilityj 
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Fig  C2.  SchesmClc  diagrams  showing  that  the  bands  of  System  B  have  been 
twisted  by  an  angle  of  ii/2  at  the  end  attached  to  the  hub. 
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X (x)  is  th6  moment  of  inertle  et  position  x  • 

Assuming  that  the  angle  twisted  is  linearly  dependent  on  x  , 


the  following  relation  can  be  obtained 
I(x)  “  cos^  X  +  sin  X 


(C-11) 


where  x  =  irx/2L^ 

Inserting  equations  (C"'9)  and  (C~ll)  into  equation  (G-10) »  one 


obtains 


twisted  •  C(I^  -  x)^/E^(I^cos2  X  +  I2  sin^  X)]  dx 

^  o 

(C-12) 

For  the  case  when  the  bands  are  not  pre-twisted,  the  value  of 
yp/Q  can  be  obtained  by  using  the  following  equation. 


^^p^^^nontwisted 


(C-13) 


The  ratio  of  the  translational  stiffness  for  the  above  two  dif¬ 
ferent  conf igutations  can  be  obtained  as  follows: 


twisted 


^^S^nontwisted  ”  ^^p'^^^nontwisted^  twisted 


Combining  equations  (C— 8)  and  (C-14),  one  can  compute 


using  the  following  equation 

/  ^^b  2  ^  ^^P^^^nontwisted 

K  .  "  12  (  — —  cos  <p  +  TTT^s - 

blP  y  ^  twisted 


(C-15) 


The  numerical  value  of  equation  (C-14)  pertient  to  the  Sandia 
System  B  design  is  approximately  equal  to  0.678.  The  angle  if  in  system 


-75- 

B  is  about  2.13“.  Using  these  values  in  equation  (C-15) ,  one  finds  that 
the  difference  of  for  both  Systems 


A  and  B  is  within  0.2S. 


APPENDIX  D 


DETEEMINAIION  OF  OUT-OF-PLANE  BAND  STIFFNESS 

To  derive  an  expression  for  the  out-ofrplane  band  stiffness,  it 
is  necessary  to  consider  a  set  of  two  bands  as  a  double  beam  with  column 
action  on  the  compression  side  of  the  beam  and  tie-bar  action  on  the  ten¬ 
sion  side.  Although  there  can  be  some  geometric  nonlinearity,  only  a 
linearized  analysis  is  presented  here.  Fig.  D1  shows  schematically,  in 
exaggerated  fashion,  the  geometric  relationships  involved. 

Denoting  the  upper  band  as  band  A  and  the  lower  one  as  band  B, 
one  can  write  the  following  expressions  for  the  respective  horizontal 
and  vertical  displacements  of  the  band-to-rim  junctions  when  the  rim  is 

tilted  through  an  angle  : 

u.,  “  R,(cos  ({)  -  1)  -  b  sin  ((I  "  -  bq) 

AL  X 


w  »  R,  sin  4  +  b(cos  -  1)  ^  Rj^ij 

iiXd  X 


(D-l) 


u  -  R.(cos  4,  -  1)  +  b  sin((i  «  biji 
BL  X 


w  “  ^  ”  b(cos  (J)  -  1)  ^ 

BL  X 


where  the  approximations  at  the  right  are  based  on  the  first-order  small- 

angle  assumptions,  sin  (J  and  1  -  cos  0  . 

The  total  axial  forces  in  the  respective  bands  are 

^A  ’  ^i  ^c  *  ■  ''aL 

^  (D-2) 

“  ^i  ^c  Vb^"BL  ^  ''bL 
^  +  A^E^(b  cos  <|)  +  Rj'  sin  (j))4>/L^ 


RIM  IN  TILTED 
POSITION 


Fig.  Dl.  Schematic  diagram  depicting  behavior  of  a  set  of  two  bands  behaving 
as  a  double  beam  with  compressive  column  action  in  the  top  beam  and 
tensile  tie-bar  action  in  the  bottom  beam.  Rotation  takes  place 


about  point  0. 
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where  P. 


=  Initial  winding  tension;  P  =  centrifugal  tension  given 


in  Appendix  B. 


The  governing  differential  equations  are 

p  .  p  -  0  (Band  A) 

dx^  ^  dx2 


(D-3) 


.4  B 

E.  ^  -  P 

dx 


.2  B 

SJL.  «  0  (Band  B) 
®  dx2 


The  boundary  conditions  are  listed  in  Table  DL 


Table  D1 


Boundary  Conditions  for  Out-of-Plane  Bending  of  Bands 


Left  End 


w  (0)  «  0 


w^(I^) 


Right  End 


cos  $  +  u^  sin  ^ 


S  (R^  cos  -  b  sin 


dw^(O) 


dw^(I^) 


„'(0)  .  0  ♦  ■  “bl  ♦>* 


B  dw®(L,  ) 

dx  dx 


in  order  to  write  the  eppropriate  solutions  for  e,uations  (D-3), 
it  is  necessary  to  know  whether  and  are  positive  (tensile)  or 

negative  (compressive).  Due  to  the  initial  tension,  it  turns  out  that 
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for  the  higher  rotational  speeds  of  major  importance,  for  small  angles, 
even  P  is  tensile  while  is  tensile  under  any  conditions.  Thus, 

the  general  solutions  may  be  .written  as 

w^(x)  *  A-  cosh  ^2  ^3  ^  ^4^ 

(D-4) 

w®(x)  »  cosh  k^x  +  sinh  k^x  +  B^  +  B^x 


where 


‘‘a  =  ‘  ^  ^2^*^ 


=  Area  moment  of  inertia  of  the  band  cross- 
^  sectional  area. about  an  in-plane  axis 


(D-5) 


Using  the  boundary  conditions  to  evaluate  A^,  .  .  *  A^  and 
•  •  *  ^  ®4*  write  the  following  results  for  the  bending  moments 

at  the  ends  of  the  individual  bands: 


-"b^2 


j2  A,-  V 
d  w  (L^) 


sinh  k.L,  A 

''aH  -  "TE - j  H  "aS.*®' 

A  b  / 

2(1  -  cosh  k^Ljj)  +  k^Ljj  sinh  k^L^ 


(D-6) 


where  R*  =  R  cos  (^  —  b  sin  iji  . 

These  same  expressions  hold  for  if  is  replaced  by  . 

Since  we  are  interested  in  the  case  where  ij>  -►  0  only,  it  is 
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accurate  enough  to  use  ^  ^  ^c  equation  (D-6) . 

In  addition  to  moments  and  associated  with  bending 

of  the  individual  bands  about  their  own  centroids,  there  is  a  moment 
carried  in  the  form  of  membrane  stresses  in  the  individual  bands.  It  is 
given  by 

-  b(Pg  -  P^)  -  2  Aj^E^jb(b  cos  ij)  +  sin  <p)<l>/^  (0-7) 

* 

The  total  moment  carried  by  the  double  beam  is 

-"“bl 

The  tilting  stiffness  of  the  double  beam  is 

^  -  M/4i  -  2  Aj^E^jb(b  cos  *  +  sin 

+  2  £^12^^ 

2  (cosh  -  1)  -  sinh 

(D-9) 

Now  it  is  necessary  to  consider  how  the  various  sets  of  bands 
contribute  to  the  stiffness  of  the  band  system.  For  a  band  located  at 
an  angle  0  from  the  normal  to  the  axis  about  which  tilting  takes  place, 
the  effective  angle  of  tilt  is  -fi  cos  0  .  Furthermore,  only  the  cos  0 

component  of  the  moment  acts  about  the  axis  of  tilt.  Referring  to  Fig. 
Cl,  one  finds  that  there  are  a  total  of  two  sets  of  bands  at  0  ,  two 
at  30*  -  0,  two  at  30*  +  0  ,  two  at  60*  -  6,  two  at  60*  +  0  ,  and  two 
at  90“  -  0  .  Thus,  the  total  out-of-plane  stiffness  related  to  band-set 

*  The  moment  generated  by  the  shear  force  where  the  band  is 
attached  to  the  rim  is  small  and  is  neglected  in  equation  (D~8) . 


sinh  k 


cosh  k^L^  - 


-) 


+  (1  -  cosh  k^L^)R* 
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stiffness  is  as  follows: 


[2  cos^-9  +  2  cos^  (30®  -  6)  2  cos^  (30®  -f  0)  +  2  cos  (60®  -  0) 


+  2  cos^  (60®  +  0)  +  2  cos^(90®  -  0)]k^  =  6  k' 


(D-10) 


It  is  interesting  to  note  that  even  up  to  0  -  32,000  rpm, 

K.  is  mainly  dominated  by  the  membrane  action,  i.e.  as  defined  in 

Dop 

equation  (D-7).  This  term  is  significantly  influenced  by  the  angle  4i  . 

In  the  Sandia  flywheel  System  A,  the  average  <5>  is  approximately  equal 
to  2.66®.  Using  equation  (D-9) .  One  obtains  23.6%  extra  stiffness  by 
incorporating  the  angle  t  as  compared  to  the  case  when  it  is  not  in- 
corporated • 

For  System  B,  the  bands  are  twisted  as  shown  in  Fig.  C2(b). 

If  a  moment. of  ,  instead  of  a  force  Q  ,  is  applied  at  one  end,  the 

moment  distribution  would  be 
M(X)  - 


The  slepe  at  the  point  of  application  of  can  be  deduced 


from  the  following  equation. 

fS  2  2  1 

twisted  “  J  Cl/Eb^^2  ^  ^1 


(D-12) 


where  2  = 

For  the  case  where  the  bands  are  not  pretwisted-,  the  value  of 

*  Another  contribution  to  the  stiffness  is  the  integrated  torsional 
stiffness  in  each  band.  However,  its  value  is  very  small  even  for  System 
B  which  has  more  torsional  stiffness  due  to  pretwisting,  and  is  neglected 


here. 
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is 


/^q)  nontwis  ted 


(1/E^^l2)dx 


(D-13) 


Combining  equations  (D-12)  and  (D-13),  one  obtains  the  following 


result  pertient  to  the  System  B  design 

j^^gj.gjj/^'^/^o^nontwisted 


(D-14) 


Thus,  the  second  term  in  equation  (D-9)  should  be  multiplied  by 

5  in  order  to  incorporate  the  pretwisted  angle  of  • 

In  System  B,  J  is  approximately  equal  to  2.13“  .  Thus, one 
obtains  approximately  19%  more  stiffness  by  incorporating  the  angle  ^ 
as  compared  to  the  case  when  it  is  not  incorporated. 
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The  method  used  here  to  calculate  the  turbine  shaft  deflection 
is  the  principle  of  virtual  work  (Casigliano's  least-work  theorem)  as  used 
for  stepped  diameter  shafts  by  Bert  [36].  Since  the  present  shaft  is 
relatively  short,  the  additional  deflection  due  to  transverse  shear  is  also 

included  here. 

From  statics,  the  bending  moment  and  transverse  shear  force 

distributions  are  found  to  be 
’m-Fz 

o  o  ■ 

(E-1) 

M<2)  “  ‘S  m  f  l_ 

(7^  '  -T—XI'  -  2)  4  1  z  1  Sb 

IS  S 

0  ^  z  1. 

S- 

where  =  S  S  ^  I-eESth  of  shaft  having  area  A^,  U  2 

overhung  distance  (between  point  of  load  and  the  closest  bearing), 

Lg  =  distance  between  bearings,  and  z  =  position  measured  from  point 

of  application  of 

The  principle  of  virtual  work  can  be  applied  to  obtain  the  following 

expressions  for  transverse  deflection  and  slope  at  the  point  where  F^ 

and  M  are  applied: 
o 

^  »  IIL  .  A  =  (E-3) 

3F^  ’  S 

where  U  is  the  sum  of  the  flexural  and  transverse  shear  strain  energies 
given*  by 


■  -S 

.  h  ^B 
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fS:B  &t(2)]^  dz 
^  “  2  E  1(2)  ^  J 

*  s 


rv(2)j  dz 

2G  K(z)A(2) 
s 


(E-4) 


A(2)  is  the  distribution  of  cross-sectional  area,  is  Young's  modulus, 

G  is  the  shear  modulus,  I(z)  is  the  distribution  of  area  moment  of 
s 

inertia,  K(z)  is  the  distribution  of  the  transverse  shear  correction 
factor,  a  is  the  distance  between  turbine  center  to  the  closest  turbine 


shaft . 


Combining  the  preceding  two  equations,  one  obtains 
mM_  V  — 

,  .  f'tB  Is  ,,,  f'tB  SS,, 

J  El  J,  K  A  G 


M  3^ 

^  aF 

O  J  _ 

'Sb 

viv. 

SF 

O 

J 

a 

K  A  G 

3 

a 

o 

S 

A 

o 

‘^B 

“  an 

‘S:b 

3H 

o 

El  ‘J*  +  j 

A 

K  A  G 

(E-5) 


(E-6) 


Inserting  the  distributions  A(2),  I(z),  M(z),  and  VCz)  approp¬ 
riate  to  the  present  general  geometry,  as  shown  in  Fig.  El,  into  equations 
(E-5),  (E-6)  and  performaning  integration,  one  obtains 


a  F  +  a, 

1  o  2  o 


F  CL-i  M 
2  o  JO 


(E-7) 

(E-8) 


where  compliances  and  are  defined  as 
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a 

1 


®1  “  ‘  4  ^  ^  ~  %  L^_ 

3  hhs  ^  ^S^BS  ^  ^S  hh^3 


•■(? 


2  2 

1  -  ^o 

^s^TS 


2  ^s^BS 


^  ^s^BS 


^  Vs 


(E-9) 


T-1 

a.  -  a  L-  -  a-  L  W 

To  calculate  the  transverse  shear  correction  coefficient,  which 
allows  for  the  nonunifom  distribution  of  transverse  shear  strain  through 
the  depth  of  the  shaft,  the  following  expression  presented  by  Cowper  [37] 
is  used 

j,  . _ 6  (1  +  v) _  (E-10) 

7  +  6v  +  4(4  +  3v)b2/(1  + 

where  B  =  inside  diameter/outside  diameter  and  v  is  Poisson's  ratio. 
E.2  Flywheel  Shaft 

The  flywheel  shaft  is  essentially  supported  as  a  cantilever  at 
its  top  end  by  rigid  attachement  to  the  turbine  shaft.  The  flywheel  shaft 
Is  quite  flexible  in  comparison  to  the  turbine  shaft  and  the  former  is 
concsutrically  located  Inside  the  latter* 

The  upper  portion  of  the  flywheel  shaft  has  a  length  of 
L  (“  10  in.)  with  an  area  moment  of  inertia  (=0.00485  in  ) .  The  lower 
portion  of  the  flywheel  shaft  contains  a  short,  necked-down  "breakaway" 

4 

region  of  length  Lj^(*1.0in)  with  area  moment  of  inertia  Ijj^(=  0.002425  in  ) 
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Uslng  the  method  of  principle  of  virtual  work  as  presented  in 
[36],  one  obtains  the  following  deflection  and  slope  for  the  hub  relative 
to  the  turbine. 

r.  I  •  O/  +  “e  \ 

h/t  4  h  5  h 

♦b/t  •  “5  +  “6  \ 

where  and  are  the  force  and  moment  applied  on  the  flywheel 

shaft  at  the  position  where  hub  located;  a^s  are  defined  as 


^  [(L^/I  )  +  (L^  -  J 

E  *“  u  u  s  u  X 


*4-3 


V  [<I-X>  <'•4  - 


‘5  ■  2  E 


(E-12) 


»6  ^  -f-  *  "-s  -  '•ui/'ti 

s 


L  a  L  +  L.  (=  11  in.) . 

S  U  X 


When  a  force  is  applied,  a  force-couple  system  .  Fj^ 

is  developed  on  the  turbine  shaft,  which  in  turn,  causes  the  turbine  shaft 
to  have  the  following  deflection  and  slope  at  the  point  where  the  flywheel 
shaft  is  attached  to  it 


“1  ^h  “2  ^s  ^h 


(E-13) 


^t  “  “2  ^h  “3  ^s  ^h 

The  net  deflection  and  slope  for  the  shaft  at  the  axial  position 
where  the  hub  is  located  (under  the  point  of  application  of  load  Fj^)  can 
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APPENDIX  F 


CALCULATION  OF  THE  VISCOUS  DAMPING  COEFFICIENT 
FOR  THE  DAMPER  IN  THE  SANDIA  SPIN-TEST  FACILITY 


In  the  Sandia  spin-test  facility,  there  is  a  Barbour  Stockwell 
damper  located  on  the  flywheel  shaft.  The  geometric  dimensions  for  this 
damper  are:  damper  radius  »  1.5  in.,  damper  axial  length  «  2  in.  Thus, 
the  total  force  acting  on  the  damper  along  the  line  of  centers  is 

r2-T 

F  -  (1.5) (2.0)  P(e)  cos  0  d9  (^"1) 


where  9  is  the  circumferential  position,  and  P(9)  is  the  total  dynamic 
pressure  acting  on  the  damper.  A  conservative  estimate  of  P(e)  is  the 
hydrodynamic  pressure,  i.e.,  the  rate  of  change  of  inlet  pressure  is 
neglected.  Following  Bansal  and  Hibner's  work  [38],  one  can  write  the 


following  expression  for  the  viscous  damping  coefficient  C^ 


^  -  [6(1.5)  y/2c] 

3r_ 


mm  mm  ry 

4  cos  9  -  c  cos  26  cos  9 
2  2 

o  (1  - e  cos  9)  (2  +  e  ) 


-6  2 

where  y  *  oil  viscosity  *  2.7  x  10  Ib-in/sec 
c  -  radial  clearance  of  damper  «  0.040  in. 
c  ■  eccentricity  ratio  (jr^l/c) 

Using  the  data  listed  above  as  input  and  integrating  equation 
(F-2)  numerically,  one  can  compute  as  a  function  of  e  .  The  results 

are  listed  below: 

£(S  |5^1/C)  0.1  0.2  0.3  0.4  0.3  0.6  0.7  0.8  0.9 

c  (Ib-sec/in)  0.109  0.114  0.123  0.139  0.165  0.21  0.295  0.497  1.3 

d 
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APPENDIX  G 


COMPUTER  PROGRAM  DOCUMENTATION  AND  LISTING 

The  numerical  computations  for  the  present  investigation  were 
carried  out  by  means  of  a  specially  written  digital  computer  program. 

The  program  was  written  in  FORTRAN  IV  language  and  was  run  on  the  IBM  Sys¬ 
tem  370/158J  at  the  University’s  Merrick  Computing  Center. 

The  program  consisted  of  two  parts.  The  first  part  was  used  to 
investigate  the  free  whirling  and  stability.  The  second  part  was  used  to 
study  the  forced  whirling.  Several  subroutines  were  used  in  addition  to 
the  main  program.  Major  variables  were  defined  by  using  comment  statements 
in  the  computer  program  itself.  The  following  is  a  list  of  the  subroutines 
and  their  functions. 

COMPLI  computing  the  compliance  coefficients 

EIGCC*  computing  the  complex  eigenvalues  of  a  general  complex  matrix 
KBAND  computing  the  total  band  stiffnesses 

LEQT2C*  solving  the  simultaneous  linear  complex  equations  or  inverting 
a  general  complex  matrix 

MULTYC  performing  multiplication  of  two  complex  matrices 


This  is  a  standard  subroutine  found  in  IMSL  Library  1.  6th  ed.. 
International  Mathematical  and  Statistical  Libraries,  Inc..  1977,  and  thus 
it  is  not  listed  herein. 
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FrttE  VIEKATICN  i  STABILITY 
IMPLICIT  REAL*a(A-H,C-2) 

H£AL«e  MR* IWtJR»MH»Ih»JH*MT .IT*JT (MD»CDt«K(32 )»IU»JO*AMPia) 
complex *16  MI  (a*8)*N(d*8)*  ALPhI 8*8 ) .K(a»8)»WA{ea)*8l(8*a) 
COMPLEX*  16  E2(e.a).C<16.16).l»{Io).Z(16.16>»F(a).A(a*8) 

COMPLEX  *16  AKbIP.  AKBOP.Al  ,  A2  .  OCMPLX  »CPM 
DATA  NOEG.NCEGO/d.ia/ 

Input  oampek  mass  ano  mass  moments  jf  inertia 

INPUT  RIM  MASS  A.NC  MASS  MOMENTS  OF  INERTIA 
INPUT  HUE  MASS  AnO  MASS  MOMENTS  OF  INERTIA 
INPUT  TURBINE  MASS  ANO  MASS  MOMENTS  OF  INERTIA 

OAT A  MO  » 10  *  JO/ 0  *0  08400 .0 .0  05900 .0*0  09500/ 

OATA  MR  *  I R  * JS/0 .04SC00  *  1  .oSoCOO  ,3.215000/ 

OATA  MH, IH. JH/0. 0235000 . J. 247000  *  0  *  0340 00/ 

OATA  MT. IT .JT/O* 02000 *0*04000 .0.067000/ 

IOAMP=0  INTERNAL  OAMPING  NOT  INCLUDED 
IOAMP=l  INTERNAL  OAMPING  INCuUCED 

I0AMP=0 

CO-VISCCUS  OAMPING  COEFFICIENT 
C0  =  0  . 

00  100  IM=l*N0EG 

00  100  IN=1 ,N0EG 
N(IM*lN)=(0.O0.0*D0) 

100  Ml  I IM.IN»  =  (0*DO*0 .DO) 

INVERT  THE  M  MATRIX 

MI  (  I  .  l)  =  l ./MR 
MI {2.2)=l ./IR 
MI (3,3)=l ./MH 
MI(4.4)=1./1M 
MI t5,5)=l ./MT 
MI  16,6) =1 ./IT 
MI (7. 7)=1 ./MO 
MI (S,8)=l ./lO 
N  (  7  *  7  )  =  OCMPLX (0. 0000. CO) 

CALL  CQMPLI lALPH, IDAMP) 

A1=ALPH(1,I) 

A2=ALPH(2,2) 

DC  900  1=1.4 

DO  900  J=I *  9 
RPM=J*lO.**i 


i 

i 


An  o 
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ftAC=KPM*3. i4l 592654/20. 

N(2.2}=JP*ftAO 
N(  4 .4  )  =  JH4KA0 
N(  e.6)  =  JT*FtA0 
N(2»6)=‘JD*RAa 

CALL  KaA^D  TO  COMPUTE  INTEGRATfeO  BAND  STIFFNESSES 

CALL  kBANDI A<clP . AKECP .RAD . lOAMP) 

ALPMI I . I )= Al+I ./AKdlP 
ALPh(  2.  =  ./AKBCP 

DO  102  IN=l,NOEG 
DC  102  JK=1 .NOEu 
IF(IK.EQ.JK)  GO  TO  101 
K(1K.JK)=(0 .00*0.00) 

GC  TC  102 

10  1  <C(  iN.uK  )=(  1  .00  *  0  .CO  ) 

102  CONTINUE 
DO  102  IA=1. NUEG 
DO  103  JA=l.iMCEG 

103  A(  lA*  jA)  =  ALFrt(  lA.  JA) 

C 

C  CALL  LEGT2C  TO  INVERT  COMPLIANCE  MATRIX 

C  LEQT2C  IS  A  ROUTINE  IN  IMSL  PACKET 

C 

C  ALL  LECT  2C(A.NDcG.NDEG.K.NDcG.NDcG.0.i*A*I  ER  } 

CALL  MwLTYCIMI  .N.ei .NOEG) 

CALL  MULTYCtMI .K. c2«NUEG) 

DC  10  KO=l.NDEG 
DO  10  CC=1 »NDEGO 
10  OI<D*LD)  =  (0  .DO.O.CO) 

00  20  KO=l.NO£G 
KCC=KC4-NO£C- 

20  D( KC.KOD)=(l .00.0.00) 

00  30  KO=l »NOEG 
KCC=KO<-NOEG 
Du  oO  LO=l  t  NOEG 
30  D(KCO.LO)=  c2(KC.LO) 

DO  40  kO=:1*NOEG 
K02=KC+NOEG 
00  40  LD=1»NDEG 
L00=L0+MJEG 

40  0(KCD.LDO)=  JlIKO.LO) 

PRI.'.T  lOOO.RPM 
PRINT  4000 
print  5000 

c  .  . 

c  CALL  EIGCC  TO  SOLVE  EIGEN  VAuUE  PH03LEM 

C  EICCC  IS  A  ROUTINE  LISTED  IN  I SML  PACKET 

C 


UUU  UVifU  ouuu  uu 


-95- 


call  EIGCC(  C.NCEGC»NCEGJ#0*i»»2.N0oG0,AK.  IER) 

DO  60  IP=i  .NOEGU 
CPM=<f(  IP  J<‘30./3  .  HI  £S2c54 
60  PRINT  2000»CPM 
900  CONTINUE 

1000  FORMAT (5X.  'KOT  SPEED ( KPM 1= • . D 15  .5 . / J 
2000  FORMAT! lOA »20 I  5 .5 } 

4000  FORMAT!  •  kE!  J(HIR  FRcl  XMI’aHIR  FRE)*) 

5000  Format !•  !Cpm)  !Cpm)*) 

STCP 

ENO 

SUEROUTINE  C3MPL I ! ALPH . 1 O AMP ) 

IMPLICIT  Cu.«PL£X*161A-0) 

CCM.PLEX*1o  ALPh(6ta) 

REAL* a  ES.  IS. ITS. I ES  tLT.LH  »LE  .LM.LO .LS » GAMAES .EIS 
RE  AL  *  a  GS  .  AO.  AI  tC  AiM  ACS.GI^.KT  .AT  .A3.Kd.LU.  lU.  IL 
C 

C  INPUT  SHAFT  aR EA  MOMENTS  CF  INERTIA 

C 

OAT A  IS.ITS.I3S/0.y04d5000,0.2375300. 0.0e6c5000/ 

C  . 

C  INPUT  SHAFT  LENGTHS  DEFINED  IN  FIG.  l.l 

C 

DATA  LT .L3 .LE .LH/l . t3D 00 . 2 .57D0  0 .  1.940  0  0 .4.06DC0/ 

INPUT  MATERIAL  PRCPEhTIES 

DATA  ES.CS/2y .SUo.l I .506/ 

DATA  GAMAES  .GAMA GS/ 0 . J  OSD  0 . O . 00 uDO/ 

INPUT  lengths  CF  AO.Al  DEFINED  IN  FIG.  El 

DATA  AO. A1/0.57D0.1  .37500/ 

INPUT  AREAS  CF  SHAFT 
INPUT  SMEAR  CCRRECTICN  FACTuH 

DATA  KT.AT,KO«AS/0.c3afcU0.1.396DJ»0.6a7';D0.0. 736300/ 
DATA  LU . 1 U.  IL/1 0 .00.0.0043500.0  .00242500/ 
Li;=LT+LE+LE 
LS-LC+LH 


. . internal  damping  is  nct  included 

.  .  . . internal  damping  13  INCLUDED 


IF! ICAMP.EG.O )  GG  TC  10 
EIS=-GAMA£S*ES 

a£s=ocmplx les.eis) 

G  I  S-=-GAM assess 
AGG^OC.MPLX!  G3.G1S) 

Go  TC  20 
10  AES=ES 
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AGS=GS 

20  ALPH1=(A1**J— A0**3) /3 . /A£  3/ 1 T  S+ ( LT*»3“A 1 ♦♦3+L  T**2*LB) /3 • 

1  /ACS/I BS+  <  A1 -AO  >/KT/AT/A0S+(LT-A1+LT**2/LB J/Ke/AS/AGS 

Al.PH2=-(  Al**2-AO**2)/2./A£S/I  TS-d-T^^a-Al^^aj/a./AES/IBS 

I-Lr*L8/3./AeS/ IGS-LT/LO/KH/AB/AGS 
AUPH3=(  Al- AO)  /  AES/ITS+{CT-A1  ) /A£5/  I  BS  <-l.a/3  ./AES/ 1  BS+l  •✓l-B/KO/AB 

1/AGS 

ALPHA'S  (  (  lU  )  ♦  (LS**3-LU**3)/'  IL)/3  ./AES 

ALPH5=(  (LU**2/Hj)  ♦•(LS*»2-LU**2)/IL)/2./AES 
ALPH6=(LU/  rL-KLS-LU  )/lL)/AES 

ALPH(l  .1  )  =  ALPH4+ALPHl*-2.  +  ALPh2»l.3+ALPHJ  *LS**2 
AUPh( 1 ,2)=ALPH5+AuPH2+AL?rt3»LS 
ALFHd  .3)  =  ALPH(  1  .  1  ) 

ALPH(1.4)=ALPH( I .2) 

ALFHd  .S}=ALPHl+ALPh2*L3 
ALPHd  .6)  =  ALPh2+ALPH3’CLS 

ALPH{  1  .  7)=L0**3/3./A£S/I  J<-ALPH1+ALPH2*l.D*(  ALPh2+ALPH3*LD) 

1 *L0+LH*C  LC»*2/2 ./ AES/ IU  +  ALPH2  FALPhB  »L3 ) 

ALPH{2.2)=ALPHfc+ALPh3 
ALPH(2.3)=ALPH( 1 .2) 

ALPH(2.4)sALPH(2.2) 

ALPH(2,S)=ALPH2 
ALFH<2.6)= ALPH3 

ALPH(2.7)= (LOAS2/2 ./AES/ IU+ALPH2+ALPH3ALC ) 

ALPH( 3. 3)  =  ALPH(  1  .  1  ) 

ALPH(3.4)=ALPH(  1 .2) 

ALPH13. 5)  =  ALPHd  ,3) 

ALFH{3»6)=ALPH( 1.6) 

ALPH(3. 7)  =  ALPHd  .7) 

ALPH(4.4)=ALPH{ 2.2) 

ALPH(4.5)=ALPH(2.S) 

ALPHtA .6)=ALPH(2.6) 

ALFH(4,7)=ALPH(2.7) 

ALPH(5.5)=ALPhl 

ALPH(5»6)=ALPh(2.S) 

ALPrilS. 7)=ALPH1 +AlPH2*L0 
AcPH{6.6)=ALPH(2.C) 

ALPH(6. 7)=ALPH2tALPH3*LO 

ALPH  C7.7)=LC*»3/3./AeS/lUFALPHl4-Al_PH2*LO-t-(  ALPh2+ALPH3*LO  )  ALO 
ALPH ( 1 • B )=  ALPH2  FL  C*  *2/2. / aES/ IJ*1 ALPH3FLD/ AES  /  lU) *LH 
ALPH(2.a)= ALPHO+LC/AES/I J 
ALPH(3.a)  =  ALPHd  .8) 

ALFH( A, 6 )= ALPH<2. fc  ) 

ALFH(5. 6)=ALPH2 
ALPHl 6. a)=ALPH3 

ALPH  {  7.  d)-=  ALPH2-)LC**‘2/  2  ./AcS/  lU 
ALPH ( b»3 )  =  ALPH3+LC/ AES/  lU 
DO  30  1=2.6 
K=l-1 

DC  30  J=1 .K 


u  u  u  u  u  w  u  u 
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30  ALPM(  1. J)-ALPH( J.I  ) 

RETURN 

£N0 

SwBHCUTINE  KaANC( AKeiP.AKdCP.RAD. IDAMP) 


XCAMP=0  .  . .  .  ,I NTfcRNAL  DAMPING  XG  NOT  INCLUDED 
IUAMP=1 . XNTERNAt.  DAMPING  IG  INCLUDED 


IMPLICIT  REAL*a( A-M.C-Z) 

CCWPLEX*lt>  CDSURT.CCeXP.OCMPLX 

COMPLEX* lo  AEBAND.D 1 .UOtUEO .FSCPR 1 t P . AK  .  AKL t AKB* AKA. AK AL. G t AM 
COMPL£X*16  hySIU.hrCOS  .AKSIP.AKoOP 
PEAL*e  LBANC 

DATA  RORl  M  .ARIM.RRI  M.ERIM/1  .3990-4*5.5:^600*8 .633D0 . 18.06/ 

2.66  IS  THE  BANC  ANGLE 

PHiiAV E=2.o8 *3 . 141592654/180  . 

TENSIO  IS  THE  EANO  INITIAL  WINDING  TENSION 
TENSIC=360.C0 

...Na=#  OF  BAND  SETS 
N0  =  12 

TO=360./Nd/2. *3. 141592654/180. 

DATA  RH*ftI/2. 2500* 7.62500/ 

LdANO=< RI-Rh)/OCC£(PhWAV£) 

C 

c  oa.ra  are  banc  width  and  thickness  respectively 

C  ROOAND.  EBANO  ARE  BAND  MATERIAL  DENSITY  AND  YOUNG'S  MODULUS 

C  GABAnD=  EANO  Less  TANGENT 

C 

DATA  SB  *T£j  .RQBANO  .EEAND/O.5D0 .0 .100 . 1  .2950-4  *  1  I  .06/ 

DATA  GABANC/3  .31 7200/ 

AaANO=aB*TO 
Al IP=Ta*da**3/l2. 

AlCP=To**34EB/ 12. 

IF! iOAMP.EG.D)  GO  TC  1 
EieAND=-GAaANO*EaANC 
AEBANDsOCMPLX ( EBANO . E 1  BAND  ) 

GC  TO  2 

1  AEBAnDsEBANC 

2  UdO=RU8ANO*KAO**2/AEBANO* ( R I ♦ *3/3-R I * *2 ♦KH/2 . ♦fiH**3/6. } 
UP=HA0**2* (KOR IM* AR IM*RH1M**2*0  TAnI TO ) ♦RC3AN0*AQAND*0CuStPHWAVE ) 

1**3* (RI**3/3.-R l**2*fiH/2 .+RH**3/6 .1/2. /LB AND) 

DR=AK  I M/RR  I  M*D  TAN  (TO  ) 

DI=A3AN0*EE AnD/2./LEAN0*DC0S( PH4AVE )**2 
UO-UP/(DR40I) 

FBCPRl=AdANC*AEBAND/LBAND*(U3-U30*OCUS( PH4AVE  )  ) *0CDS ( PH* AVE 1 
f aC=HL8ANO*ABAND*KA04*2*C  2  .♦«  !♦ *3-3  . ♦R I  * ♦2*RH ♦KM**3 ) /6 . /I R I-RH ) 
P=TENSIC+FeCRFECPRl 
AK-=CDSORT(P/AEEANC/Ar  IP) 
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AKL=AK#UclANC 

HYSlN=(CDexP(AKL>-CC£XH(-AKL)  )/2. 
hVCOS  =  (CDEXF( AKL ) ♦CD£XP(-AKLJ  )/2. 

AKa=AK*F*HYSlN/(2.*< 1,-HYCCS) ♦AKC»HYSIN) 

AK8IP=12 ( ASANC*  AEfi ANO/LeANO*DCCS (PmW A VE } ♦♦2+AK8 ) 
AKA-=COSQHT{F/A£aANO/A  ICP) 

akai.=aka*lband 

HYSiN=<COEXP( AKAL)-CC£XP(-AXAL) )/2* 

HYCOS-=(  COEXP(  AKAL  )+CO£XP (  -AKAL)  i/2» 

RP=RI ♦OCOS (PHto AVe )- 1 .£*03  IN (Phw AVE) 

U=CHYC0S-HYEl4N/AKAL>*LbANi:+(  I  .-HYCC3JARP 
A»»=A£3AN0*  a  IDP*  AK  A**2AG/  (  2  .♦(  HYCu5— I  •  )- AKAL  AH  Y  SIN) 

AKa=2  .♦  AdANC*  AEDANO+I  •£♦  (1  ,b*2CCSI?Hi»AV£)+RI*CSIN(PHWAVE)  )/LBANO 
14-2. AAM 
AKeCP=6.*AKE 
KETUKN 
ENC 

StidPCUriNE  kULTYCtA»B»CtN» 

C0»PLeX*16  A(N«N) ♦a<N.N)»C(N*N) 

OQ  B  1=1  xN 
DC  B  U=1*N 
Cl  l»J)=((). 1)0*0.00) 

00  5  K=1 .N 

5  C(  I  .  J)=CI  1  .  J)  t-At  I  .K)*e  (K.J) 

RETURN 

bNC 


nr>r>  nnn  nf>nf>n  cnnnno  nnr»A  nr>n 
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FORCEO  VIBRATION 

IMPLICIT  REAL*a(A-H.C-Z) 

|iEAL*8  MR,  Ifi.JR.Mh.Ih.JH.MT.IT,  JT,MO«CO,»K(.32).IO.JO.AMP{a) 

COMPLEX*  1 1»  ALPHI  a  »a)  .K(3»a)«F<a)  »«fA(8a)  »A(6.8) 

COMPLEX *16  OCMPLX*AKB IP* AKBCP  tA 1 • A2*CPM 

IOAMP=il  INTERNAL  DAMPING  NOT  INCLUOED 

lOAMPsi  Internal  damping  included 

IDAMP=l 

DATA  N0£G*NCeG0/3* 16/ 

INPUT  DAMPER  MASS  AND  MASS  MOMENTS  OF  INERTIA 

Input  rim  mass  and  mass  moments  of  inertia 
INPUT  hua  mass  and  mass  moments  cf  inertia 
input  turbine  mass  and  mass  moments  of  inertia 

DATA  mO*ID*JO/0.0034C0.0.0C59D0. 0.009500/ 
data  HR.IR.JR/0.043D00.1 .636000 . 3 .21 5000/ 
data  MH.IH.Jri/3. 0235000.0.247000 .0.354030/ 

DATA  MT .IT . JT/0 .02000 » 0 .04000 .0.067000/ 

I  c  A  3  E  =  1  initial  unbalance 

I  C  A  S  E  =  2  initial  TILT 

I  C  A  S  E  =  3  COMBINATILN 

ICASc=l 

DO  900  ICO^I.IO 

C0=VISCCU3  DAMPING  COEFFICIENT 

CCsICC*0.l 
PRINT  7000.  CO 
IF(lCASE-2>  S.6.7 

5  PRINT  4000 
GC  TC  6 

6  PRINT  5000 
CO.TC  8 

7  PRINT  6000 
3  CONTINUE 

call  COMPLI  TC  COMPOTE  COMPLIANCE  COEFFICIENTS 

CALL  CCMPLI (AL PH, lOAMP) 

AI=ALPM( 1. I  ) 

A2=ALPHI2»2) 

00  900  1=1.4 
DO  900  J=l .9 
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Rpy=j*i 0 .♦* I 
F(A0=RPM*3.  141  592654/30* 

RA0=RPM*3. 141592654/30. 

CALL  KBAND  TO  COMPUTE  INTEGRATED  BAND  STIFFNESSES 

CALL  KBANOI AK8 IP  »  AKECP.R AO. I  CAMP } 

AL?H( 1»1)=A1+1./AK31P 
ALPH(2.2)=  A2<-1  ./AKBCP 
00  102  IK-l.NDEG 
DO  132  JK=1,N0EG 
IF(IK.EQ.JK)  GO  TO  101 
K(  IK. JlO  =  (  0.00.0.00) 

CO  TO  102 

101  K( IK. JK )  =  {  1  .00 .0  .00  ) 

102  continue 

00  103  lAsl.NDEC 
DO  1 03  3A=1.N0EG 

103  ACIA.JAI^AL.  UIAtJA) 

CALL  LEQT2L  TO  INVERT  COMPLIANCE  MATRIX 
LE0T2C  ISA  ROUTINE  In  IMSL  PACKET 

CALL  LEQT2C (A .NOEG. NDEG.K .NOEG. NOEG .0 .W A. lER ) 

DO  104  1F=1»NDEG 

104  F( IF)  =  ( 0.00.0  .CO) 
lF(ICA5E-2)  105.106.107 

105  r  (  1  )  =  MR4HA0»42 
GC  TO  lOa 

106  F(2)=( IR-JR)*RA04*2 
GO  TO  106 

107  F( 1 )=  MR*RAD**2  *0.028 

FI2)=  ( IR- JR ) *RAD**2*0 .0023 

lOa  K(1.1)=KU  .n-MR*RAC**2 

Kl2»2)=K(2.2)-< IK+JK )*RAO**2 
K(3»3>-K(3 .3 )-MH*RAD**2 
K(4«4)=K(4.4)-( IH+JH) *«AD**2 
K(5.5)=K(5 .5)-MT*RAC**2 
K(6.e)=K(6.e)-(lT+JT )*RAD**2 

KI 7,7)=K( 7.7)-M0*KAC**2-0CMPLX( 3 .0000 .CO )*K AD 
K(8.8)=K<a.tt)  -( IC*JO) ♦RAC**2 

CALL  LECT2C  TO  SOLVE  LINEAR  SIMULTANEOUS  ECS 

CALL  lEGT2C(K.N0£G.NCEG.F. 1 .NUEG.O.aA. XER) 

DC  201  IP= 1 .NDEG 
201  AMP(IP)=COAeS(F(IP))  . 

PRINT  3000  .fiPM. (AMP(L ) .L=l »nO£G) 

900  CONTINUE 
3000  FORMAT { 1  OX .901 2 .4 ) 


uou  uuu  uou  uuuu 


-101- 


4000  .^CKMAIC  RCT  SPEED  ASS  RR/ER  ABS  PR/ER  A8S  RH/ER 

1  ABS  PH/cR  ABS  RT/ER  ABS  PT/ER  ASS  RO/EF  A8S 

5000  FORM AT (•  RCT  SPEED  ABSRR/PCR  ABSPR/PCR  ABSRH/PCR 

1  ABSPH/POR  ABSRT/PCR  ABSPT/PUR  ABSRD/PCK  ABSPO/POR  ’  »/ ) 

6003  FCRMATC  RCT  SPEED  ABSRR  A8SPR 

1  A8SPH  A8SRT  ABSPT  ASSRC  A8SP0  t/) 

7000  F0RMAT(5X, *00=* *F5.2./) 

STCP 

END 

SUaPOUTINE  CQMPLi  (Ai.PH^iDA^'P) 
implicit  CCMPL£X416(A-D) 

C0MPLEX416  ALPhie.a) 

KEAL*d  ES. IS. lT5.IfcS.LT.ue.LE.CH.L0.LS.GAMAES.EI3 
hgAL+a  CS.AO.AI.GAMACS.uIS.KT .at. ab.kb . lu. iu. il 

"  INPOT  SHAFT  AREA  MOMENTS  OF  INERTIA 

DATA  IS. ITS. 18 5/ 3. 0 04 6 SO 30.3. 2oT5D 03. 0.3 86 65 000/ 

INPUT  SHAFT  LENGTHS  DEFINED  IN  FIG.  I.l 

DATA  uT .L3 .LE.LHy  I .63000 . J .37000.  I . >400 C . 4 . 060 CC/ 

INPUT  MATERIAL  PROPERTIES 

OATa  cS.GAMAeS/29. 5036 .0.005003/ 

OATA  u5 .GAMAGS/1 I  .506. G.0C5D0/ 

input  LENGTHS  OF  A3.A1  DEFINED  IN  FIG.  El 

DATA  AO.AI /O. 5700 .1 .37500/ 


INPUT  AREAS  CF  SHAFT 

input  shear  CORRECTION  FACTOR 

OATA  K.T  .AT  .K8.  AO/ 0  .63  5600.  I  .29600.0  .587  500.0. 736300/ 

DATA  t.U  .  IU  .  IL/1 '3 . 00 .0  •  0046503  .0 . 00242503/ 

LD=LT4LB+LE 

LS=LDHLH 

. . internal  oamping  is  not  InCuUDED 

I  . .  •  1DAMP=1  . . .  .  . I NTERNAL  damping  IS  INCLUDED 
iF(  10AMP.EG.3 »  OC  TO  10 
EIS=--GAMA£5«^ES 
AES=OCM?LX (ES . El S  > 

GIo=-GAMAG£4GS 
AG3=CCMPCX (G5. GIS I 
GO  TO  20 
10  AE3=ES 
A  O  S  ^  0  ^ 

60  ALPril=( A1443-AC**3)/3./AES/lTS4(LT*v3-Al**3+Lr**24LB)/3. 
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I  /AES/IBS+<A1-A0)  /KT/AT/AGS«-{  UT-Al  4-LT«*2/UQ)/KE/AB/A*;S 
AUPH2=  — (  A1  **2  — A0<=*2)/2»/AES/I  TS— (LT  ♦♦2— A  l**Z)  /  2  •/AES/  IBS 
l-LT*l-B/3./A£S/l  JS-LT/ua/KQ/AS/AGS 
ALPH3  =  <  Al-AO) /AES/I TS+ (LT-Al J /AE3/ I BS +Ua/3 ./AES/ 133  + 1 ./LB/KB/AB 

l/AGS 

ALFHA=( (L0++3/IU) ♦ (LS** 3-LU**3 )/ IL ) /3 ./AES 
Al.PHt>=(  (LJ**2/1U)  +  (LS*»2-LG  +  A2)/IL)/2./AES 
ALPH6-=(  LU/  I  U+  (LS-CO  J  /  ID/AES 

ALPH(  I  •  1  )  =  ALPliA  +  ALPhl  +2.*ALPh2  +  US*AUPHJ*LS  +  *2 
ALPh( 1 «  2 )-ALPH5  +  ALPh2  +  ALPH3*LS 
AL.PH  (1.3)  =  ALPH  (  1.1) 

ALPh(  1.4)  =  ALPH( 1 .2) 

AL?H( 1.5)=ALPHl+ALPh2*LS 
AuPhC  1 .6)=ALPH2«-ALPh3*L3 

ALPhi  1 .7)  =  LB»*5/3  ./ AES/ 1  U  +  ALPH  1  ♦-ALPH2*L0  +  (  AUPh2+ALPH3*LD) 
l+LD+UM*(LD**2/2./AES/ 10+ALPH2+ALPH3+UD) 

ALPH( 2.2)= ALPH6 RALPHS 
ALPh(2.3)=ALPH{ I .2) 

ALPM( 2.4)=ALPrt( 2.2) 

ALPHl 2. 5 )=ALPH2 
ALPH(2.6)=AUPH3 

ALPH(2. 7)= (LD+*2/2./Ae5/I0+AUPH2+ALPH3»L0) 

ACPH(3.3 )=ALPH( 1 . 1 ) 

AUPH(3.4)=ALPh( 1 .2) 

ALPH(3.5>=ALPH( 1 ♦ 5) 

ACPH(3.6)=ALPh{1 ,6) 

AUPH( 3. 7)=ALPH( 1 . 7) 

ALPH ( 4.4 )= ALPH (3.2) 

ALPH( 4.5)=ALPH(2 .5) 

AUPH(4.6)=ALPH(2.6) 

ALPH (4.7)= ALPH (2. 7) 

ALPH( 5.5)=ALPH1 
ALPH(S.6)=ALPH(2.5) 

A4-PH(5.  7)=ALPH1+ALPH2  +  L0 
ALPH( 6.6 )=ALPh (2.6) 

AUPH(6.7)=ALPH2+ALPH3=LU 

AUPH(  7.7)  =CU*=  3/3  ./AES/ 1  iJ+ ALPHl  +ALPH2  *LD  +  (  ALPh2+ALPH3*L0  )  *L  C 
ALPH  (  1  .8  )  =  ALPH2  +l.C*»2/2  ./AES/IU  +  (  ALPH3+LD/AcS/IU  >  ♦LH 

ALPH(2.a)=ALPH3+LC/AES/IU 

ALPh(3. 8)=ALPH( 1 .8) 

ALPH (4.8)= ALPH (2. a) 

ALPH( 6. 6)=ALPH2 
ALPH (6. 3 )= ALPH3 

AL?H( 7.3)=ALPH2+H;+*2/2./A6S/ lU 
ALPH( a.6)=ALPH3+LC/AES/IU 

OG  30  1=2.8 
<=1-1 

DO  30  J=1.K 

30  ALPH( I . J)  =  ALPH( J.  I  ) 

RETURN 
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SUBROUT  INE  KBAMC(  AKti  IP,  AKOCP.HAO.  I  DAMP) 


IDAMPsQ . INTERNAL  DAMPING  IS  NOT  INCLUDED 

10AMP=l .....INTERNAL  DAMPING  IS  INCLUDED 


IMPLICIT  REAL*8( A-H, C-Z) 

CO^«PLEX*lt>  CDSQRT.CDEXP.UCMPLX 

COMPLEX *16  AEBAnO.DI , UO » U BO . FBCPR I ,P , AK , AKL . A KB , AKA . AX AL .G , AM 
COMPLEX+lb  HTSIN.hYCCS .AKBIP, AKBOP 
REAL*e  LBAND 

2.66  IS  THE  BAND  ANGLE 

PHUAVE=2.66^3. 141592654/ laO, 

DATA  fiORlM,ARI.v,RRlM,£RIM/l  .3990-4,5.5960  0,6.63300. 18.0  6/ 

TENSIO  IS  THE  BANC  INIT IAL-  H INOINu  TENSION 

TENSIC=36D.OO 
C... . . . .  .Ne=«  OF  BAND  SETS 

Nti  =  12 

TO=360./NB/2. 43. 141592654/180. 

DATA  RH,RI/2. 2500. 7, 62500/ 

LeAN0=(  R I-RF) /DCC 5 (PH A A VE) 

C 

C  BB.T8  ARE  BAND  A  ICTH  AND  THICKNESS  RESPECTIVELY 

C  RCUAND,  EdANO  ARE  BAND  MATERIAL  CEnS  ITY  AND  YOUNG'S  MODULUS 

C 

DATA  EB.TB  .ROBAND .EEANC/0 .500 ,0 . IDO , 1 .2950-4 , 1  I .06/ 

C 

C  GADANO=  BAND  LCSS  TANGENT 

C 

DATA  GABANO/O .017200/ 

ABAND=ba*TB 
AI IP=TB*a34*3/12. 

A £CP=r 8*43  4  88/ 12 . 

IFdDAMP.EC.O)  GC  TO  1 
EI8AN0=-GA8AND4£eANC 
AEBANO=OCMPLX (EBANO.EIBANO ) 

GO  TO  2 

1  AEBANO=EBANO 

2  UdO=RCaANL)*fiAu**2/A£EAND4  (R  I4*3/3-R  I442*BH/2. +fiH443/6.  ) 
CP=«A04*2* ( kaRIM4  ARl M4RkIM4  4  24  0rAN t TO  ) HRCB AND 4 AB AND«0C0S ( P h 4 AV E ) 

1 *4341 R I  443/3. -RI4424KH/2.+RH4 43/6. )/2./LeAN0) 

0R  =  AKIM4£R  l.y/RRIM4OTAN(T0  ) 

OI=A8AN04£e AND/2 ./LB AN040C0S < PH * A VE  >44  2 
UO=UP/(OR+C1) 

FBCPRl  =  AdANC#AEQAN0/LBAN04( U0-UB040C03( PhfcAVE)  )  40C0S ( PHX AVE  ) 
FaC=t<CBAN04A8AN0  4RAD442  4(  2.  4R  14  43- 3 . 4k  I  4»24RH +RH4*3  ) /6  . /(  R  l-RH  ) 
P=T£NSI0+F6C+FeCPfiI 
AK=CDSQRT { P/AEEANC/A 1  IP ) 
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AKL=AK«:LeANO 

hVSlN=(  COEXP(  AK,L)-CDEXP(-AKLJ  )/ 2  • 

MYCOS-(CDeXP( AKLJ+COexPC-AKE) )/2. 

AKB=AK*P*HYaiN/C2.*( 1 HYCOS) 4AKL*HYiIN) 

AKBIP=12  .*  t  AB AN0»AeeAND/L0AN0*0C03(PHWAVE) •♦2+AKa J 
AKA-COSQRKP/AEcANC/AICP  ) 

AKAL-AKA*L6AN0 

HY£lN-=(COEXPi  AKAL)-CCEXP(-AKAl.)  )/2. 

HYCOS=(COEXP( AKAL )4COEXP{-AKAU) )/2. 

RP=R£*DCOS (PHW Ave )- I .5*05 IN(PbwAVfc') 

C=(HYCOS-MYSir4/AKAL)4LoANC*(  1  .-HyCCS}»RP 

AM=  AEbANO^A  10P*AKA**2»C»/(  2  .♦  {  HYCC5-  I  •  )-AKAL*H  VSIN) 

AKe=2.*AaANCCAfceAN0*  I  «S*(  I  ,‘3*CCJS(Phl*AV£  ) +H I  *  0  £  IN  ( PHKi  A  V£  )  )/uaAND 
142. ♦AM 
Ai<aCP=6.4AKc 
RETURN 
END 


